4 · Sampling and Elementary Statistical Inference
This topic follows CB, Chapter 5, Sections 5.1–5.4. Section 5.5 (convergence concepts) is Topic 5.
Random samples
A data set is a collection of observed numbers. We treat these numbers as realizations of random variables, and the simplest model for how they came about is the random sample.
Definition 5.1.1 (iid sample). The random variables \(X_1, \ldots, X_n\) are called a random sample of size \(n\) from the population \(f(x)\) if \(X_1, \ldots, X_n\) are mutually independent and the marginal pdf or pmf of each \(X_i\) is the same function \(f(x)\). The common term for such a sample is iid: independent and identically distributed.
By independence, the joint pdf or pmf of the sample is the product of the marginals: \[ f(x_1, \ldots, x_n) = f(x_1) f(x_2) \cdots f(x_n) = \prod_{i=1}^{n} f(x_i). \]
We will usually assume that the population belongs to a parametric family \(f(x \mid \theta)\), and then \[ f(x_1, \ldots, x_n \mid \theta) = \prod_{i=1}^{n} f(x_i \mid \theta). \]
Example 5.1.2. Let \(X_1, \ldots, X_n\) be a random sample from an exponential population with mean \(\beta\), for example the time until failure for \(n\) identical units. The joint pdf of the sample is \[ f(x_1, \ldots, x_n \mid \beta) = \prod_{i=1}^{n} \frac{1}{\beta} e^{-x_i/\beta} = \frac{1}{\beta^n} e^{-(x_1 + \cdots + x_n)/\beta}, \qquad x_1, \ldots, x_n > 0. \]
The probability that all \(n\) units last more than 2 time units is \[ \begin{aligned} P(X_1 > 2, \ldots, X_n > 2) &= \int_2^{\infty}\!\!\cdots\!\int_2^{\infty} \prod_{i=1}^{n} \frac{1}{\beta} e^{-x_i/\beta}\,dx_1 \cdots dx_n \\ &= \prod_{i=1}^{n} \int_2^{\infty} \frac{1}{\beta} e^{-x_i/\beta}\,dx_i = \left(e^{-2/\beta}\right)^n = e^{-2n/\beta}. \end{aligned} \] The multiple integral splits into a product of one-dimensional integrals because the integrand is a product and the region is a product of intervals. If \(\beta\) is large relative to \(n\), the probability is close to one.
Statistics and unbiasedness
The basic problem of parametric, frequentist statistics is this: we assume that the population belongs to a parametric family, but the value of the parameter is unknown.
\(X\) is the number of successes in \(k\) Bernoulli\((p)\) trials. Then \(X \sim \text{Bin}(k,p)\), but \(p\) may be unknown.
\(X\) is the time to the first event in a Poisson process with rate \(\lambda\). Then \(X\) is exponential with rate \(\lambda\) (mean \(1/\lambda\)), but \(\lambda\) is unknown.
The question is whether the observed sample \(X_1, \ldots, X_n\) can be used to find a plausible value of the unknown parameter. One answer is to choose the value of the parameter that maximizes the probability of the observed sample. This is maximum likelihood estimation, the subject of Topic 6.
Statistic. Let \(X_1, \ldots, X_n\) be a sample. Any function \(T = T(X_1, \ldots, X_n)\) of the sample is called a statistic.
Unbiasedness. Let \(X_1, \ldots, X_n\) be a sample from \(f(x \mid \theta)\). A statistic \(T = T(X_1, \ldots, X_n)\) is an unbiased estimator of \(\theta\) if \(E(T) = \theta\).
The sample mean. Let \(X_1, \ldots, X_n\) be a random sample from a population with \(E(X) = \mu\), and let \(\bar{X} = \frac{1}{n}\sum_{i=1}^{n} X_i\). By linearity of expectation, \[ E(\bar{X}) = E\left(\frac{1}{n}\sum_{i=1}^{n} X_i\right) = \frac{1}{n}\sum_{i=1}^{n} E(X_i) = \frac{1}{n} \cdot n\mu = \mu, \] so \(\bar{X}\) is an unbiased estimator of \(\mu\).
“The mean” is a concept for both the sample (the average) and the population (the expectation \(E(X)\)). The average is an observable property of the sample, and we use it to learn about an unobservable property of the population. Working out such correspondences precisely is statistical inference.
Sample variance and sampling distributions
Definition 5.2.3. The sample variance is the statistic \[ S^2 = \frac{1}{n-1} \sum_{i=1}^{n} (X_i - \bar{X})^2, \] and the sample standard deviation is \(S = \sqrt{S^2}\).
\(\bar{X}\) carries information about \(\mu\), and \(S^2\) carries information about \(\sigma^2\). The divisor \(n-1\) rather than \(n\) is explained in Why \(n-1\) below.
Definition 5.2.1 (second half). The probability distribution of a statistic \(T = T(X_1, \ldots, X_n)\) is called the sampling distribution of \(T\).
A statistic is a function of random variables, so it is itself a random variable, and it has a distribution. The figure shows the sampling distribution of \(\bar{X}\) for samples from the exponential distribution with mean 1, which has \(\mu = \sigma^2 = 1\) and is strongly skewed. Each panel is a histogram of \(\bar{X}\) computed on 20,000 independent samples of size \(n\), drawn on the same vertical scale. The dashed curve is the \(N(1, 1/n)\) density.
Two things change as \(n\) grows.
- The spread shrinks. The standard deviation of \(\bar{X}\) is \(\sigma/\sqrt{n}\) (Theorem 5.2.6 below): \(0.71\), \(0.32\) and \(0.14\) in the three panels.
- The shape approaches the normal, even though the population is far from normal. This is the central limit theorem, Topic 5.
The rest of this topic works out sampling distributions exactly: first the mean and variance of \(\bar{X}\) and \(S^2\) for any population, then the full distributions of \(\bar{X}\), \(S^2\) and related statistics when the population is normal.
Sums of random variables from a random sample
Theorem 5.2.4. Let \(x_1, \ldots, x_n\) be any numbers, and \(\bar{x} = (x_1 + \cdots + x_n)/n\). Then
- \(\displaystyle \min_a \sum_{i=1}^{n}(x_i - a)^2 = \sum_{i=1}^{n}(x_i - \bar{x})^2\),
- \(\displaystyle (n-1)s^2 = \sum_{i=1}^{n}(x_i - \bar{x})^2 = \sum_{i=1}^{n} x_i^2 - n\bar{x}^2\).
Part (a) says that \(\bar{x}\) is the least-squares fit of a single constant to the data.
Both parts rely on the fact that deviations from the mean sum to zero: \[ \sum_{i=1}^{n}(x_i - \bar{x}) = \sum_{i=1}^{n} x_i - n\bar{x} = n\bar{x} - n\bar{x} = 0. \]
Proof. (a) Write \(x_i - a = (x_i - \bar{x}) + (\bar{x} - a)\) and expand the square: \[ \begin{aligned} \sum_{i=1}^{n}(x_i - a)^2 &= \sum_{i=1}^{n}(x_i - \bar{x})^2 + 2(\bar{x} - a)\underbrace{\sum_{i=1}^{n}(x_i - \bar{x})}_{= 0} + n(\bar{x} - a)^2 \\ &= \sum_{i=1}^{n}(x_i - \bar{x})^2 + n(\bar{x} - a)^2. \end{aligned} \] The first term does not depend on \(a\), and the second is non-negative and equal to zero exactly when \(a = \bar{x}\).
- Expand the square and use \(\sum_i x_i = n\bar{x}\): \[ \begin{aligned} \sum_{i=1}^{n}(x_i - \bar{x})^2 &= \sum_{i=1}^{n} x_i^2 - 2\bar{x}\sum_{i=1}^{n} x_i + n\bar{x}^2 \\ &= \sum_{i=1}^{n} x_i^2 - 2n\bar{x}^2 + n\bar{x}^2 = \sum_{i=1}^{n} x_i^2 - n\bar{x}^2. \qquad \blacksquare \end{aligned} \]
Lemma 5.2.5. Let \(X_1, \ldots, X_n\) be a random sample from a population and let \(g(x)\) be a function such that \(E g(X_i)\) and \(\operatorname{Var}(g(X_i))\) exist. Then \[ E\left(\sum_{i=1}^{n} g(X_i)\right) = n\, E g(X_1) \qquad\text{and}\qquad \operatorname{Var}\left(\sum_{i=1}^{n} g(X_i)\right) = n \operatorname{Var} g(X_1). \]
Proof. The first result is linearity of expectation together with \(E g(X_i) = E g(X_1)\) for every \(i\), since the \(X_i\) have the same distribution. Independence is not used.
For the second, write \(\gamma = E g(X_1)\). Then \[ \begin{aligned} \operatorname{Var}\left(\sum_{i=1}^{n} g(X_i)\right) &= E\left(\sum_{i=1}^{n} \bigl(g(X_i) - \gamma\bigr)\right)^2 \\ &= \sum_{i=1}^{n} E\bigl(g(X_i) - \gamma\bigr)^2 + \sum_{i \neq j} E\Bigl[\bigl(g(X_i) - \gamma\bigr)\bigl(g(X_j) - \gamma\bigr)\Bigr]. \end{aligned} \] Each term in the first sum is \(\operatorname{Var} g(X_1)\). For \(i \neq j\), \(X_i\) and \(X_j\) are independent, so the expectation of the product is the product of the expectations (Theorem 4.2.10), and each factor has expectation zero. Here independence is used. \(\qquad \blacksquare\)
Theorem 5.2.6. Let \(X_1, \ldots, X_n\) be a random sample from a population with mean \(\mu\) and variance \(\sigma^2 < \infty\). Then
- \(E(\bar{X}) = \mu\),
- \(\operatorname{Var}(\bar{X}) = \sigma^2/n\),
- \(E(S^2) = \sigma^2\).
Proof. Part (a) was shown above. For (b), a constant comes out of the variance squared, and Lemma 5.2.5 with \(g(x) = x\) gives \[ \operatorname{Var}(\bar{X}) = \operatorname{Var}\left(\frac{1}{n}\sum_{i=1}^{n} X_i\right) = \frac{1}{n^2}\operatorname{Var}\left(\sum_{i=1}^{n} X_i\right) = \frac{1}{n^2} \cdot n \sigma^2 = \frac{\sigma^2}{n}. \]
For (c), use Theorem 5.2.4(b), together with \(E X^2 = \sigma^2 + \mu^2\) and \(E\bar{X}^2 = \operatorname{Var}(\bar{X}) + (E\bar{X})^2 = \sigma^2/n + \mu^2\): \[ \begin{aligned} E(S^2) &= E\left(\frac{1}{n-1}\left[\sum_{i=1}^{n} X_i^2 - n\bar{X}^2\right]\right) = \frac{1}{n-1}\left(n E X^2 - n E \bar{X}^2\right) \\ &= \frac{1}{n-1}\left(n(\sigma^2 + \mu^2) - n\left(\frac{\sigma^2}{n} + \mu^2\right)\right) = \frac{1}{n-1}(n-1)\sigma^2 = \sigma^2. \qquad \blacksquare \end{aligned} \]
So \(S^2\) is an unbiased estimator of \(\sigma^2\). With the divisor \(n\) instead of \(n-1\), the expectation would be \(\frac{n-1}{n}\sigma^2\).
Unbiasedness is not enough
Unbiasedness alone does not single out an estimator. For a random sample with mean \(\mu\) and variance \(\sigma^2\), the three statistics \[ T_1 = \bar{X}, \qquad T_2 = \frac{X_1 + X_2}{2}, \qquad T_3 = X_1 \] are all unbiased for \(\mu\), with variances \(\sigma^2/n\), \(\sigma^2/2\) and \(\sigma^2\). \(T_2\) and \(T_3\) ignore most of the data, and their larger variances show it. A criterion that takes both location and spread into account is the mean squared error.
Mean squared error (CB §7.3.1). The mean squared error of an estimator \(T\) of \(\theta\) is \[ \text{MSE}(T) = E\bigl((T - \theta)^2\bigr). \] The difference \(E(T) - \theta\) is the bias of \(T\).
Bias–variance decomposition. \(\text{MSE}(T) = \operatorname{Var}(T) + \bigl(E(T) - \theta\bigr)^2.\)
Proof. Add and subtract \(E(T)\): \[ \begin{aligned} E\bigl((T - \theta)^2\bigr) &= E\Bigl(\bigl(T - E(T)\bigr) + \bigl(E(T) - \theta\bigr)\Bigr)^2 \\ &= E\bigl(T - E(T)\bigr)^2 + 2\bigl(E(T) - \theta\bigr)\underbrace{E\bigl(T - E(T)\bigr)}_{= 0} + \bigl(E(T) - \theta\bigr)^2 \\ &= \operatorname{Var}(T) + \bigl(E(T) - \theta\bigr)^2. \qquad \blacksquare \end{aligned} \]
For an unbiased estimator the MSE is the variance, so among \(T_1\), \(T_2\) and \(T_3\) the sample mean has the smallest MSE. A biased estimator can still have a smaller MSE than an unbiased one, if the bias is paid for by a reduction in variance. The section The MSE of \(S^2\) and \(\hat\sigma^2\) gives an example.
The mgf of the sample mean
Theorem 5.2.7. Let \(X_1, \ldots, X_n\) be a random sample from a population with mgf \(M_X(t)\). Then the mgf of the sample mean is \[ M_{\bar{X}}(t) = \left(M_X(t/n)\right)^n. \]
Proof. \[ \begin{aligned} M_{\bar{X}}(t) = E e^{t\bar{X}} &= E e^{t(X_1 + \cdots + X_n)/n} \\ &= E\left[e^{(t/n)X_1} e^{(t/n)X_2} \cdots e^{(t/n)X_n}\right] \\ &= E e^{(t/n)X_1} \cdots E e^{(t/n)X_n} \qquad \text{(by independence)} \\ &= \left(M_X(t/n)\right)^n \qquad \text{(identical distributions).} \qquad \blacksquare \end{aligned} \]
When the result is a recognizable mgf, this identifies the sampling distribution of \(\bar{X}\). For the normal population, \(M_X(t) = \exp(\mu t + \sigma^2 t^2/2)\), and \[ M_{\bar{X}}(t) = \left(\exp\left(\frac{\mu t}{n} + \frac{\sigma^2 t^2}{2n^2}\right)\right)^n = \exp\left(\mu t + \frac{(\sigma^2/n)\, t^2}{2}\right), \] which is the mgf of \(N(\mu, \sigma^2/n)\).
Sampling from the normal distribution
Under normal sampling, the sampling distributions of the statistics that are used most in inference can be derived exactly. By the central limit theorem (Topic 5), the results remain approximately correct in much more general situations. There are three tasks:
- the distribution of \((n-1)S^2/\sigma^2\),
- the distribution of the \(t\)-statistic \((\bar{X} - \mu)/(S/\sqrt{n})\),
- the distribution of the ratio of two sample variances.
The arguments are constructive and long, but not advanced. They use the transformation techniques from Topic 3.
The \(\chi^2\) distribution
Recall that \(\chi^2_p\), the chi-squared distribution with \(p\) degrees of freedom, is the Gamma\((p/2, 2)\) distribution, with pdf \[ f(x) = \frac{1}{\Gamma(p/2)\, 2^{p/2}}\, x^{p/2 - 1} e^{-x/2}, \qquad x > 0, \] mean \(p\), variance \(2p\) and mgf \((1 - 2t)^{-p/2}\) for \(t < 1/2\).
Lemma 5.3.2 (facts about \(\chi^2\) random variables).
- If \(Z \sim N(0,1)\), then \(Z^2 \sim \chi^2_1\).
- If \(X_1, \ldots, X_n\) are independent and \(X_i \sim \chi^2_{p_i}\), then \(X_1 + \cdots + X_n \sim \chi^2_{p_1 + \cdots + p_n}\).
Proof. (a) From Examples 2.1.7 and 2.1.9 in CB, if \(Y = X^2\) then \(f_Y(y) = \frac{1}{2\sqrt{y}}\bigl(f_X(\sqrt{y}) + f_X(-\sqrt{y})\bigr)\) for \(y > 0\). With \(f_X\) the standard normal pdf, \[ f_Y(y) = \frac{1}{2\sqrt{y}} \cdot \frac{1}{\sqrt{2\pi}} \left(e^{-y/2} + e^{-y/2}\right) = \frac{1}{\sqrt{2\pi}}\, y^{-1/2} e^{-y/2}, \] which is the \(\chi^2_1\) pdf, since \(\Gamma(1/2) = \sqrt{\pi}\).
- By independence, the mgf of the sum is the product of the mgfs (Theorem 4.6.7): \[ M_{X_1 + \cdots + X_n}(t) = \prod_{i=1}^{n} (1 - 2t)^{-p_i/2} = (1 - 2t)^{-(p_1 + \cdots + p_n)/2}, \] which is the mgf of \(\chi^2_{p_1 + \cdots + p_n}\). \(\qquad \blacksquare\)
Together, (a) and (b) give: if \(Z_1, \ldots, Z_p\) are iid \(N(0,1)\), then \(Z_1^2 + \cdots + Z_p^2 \sim \chi^2_p\).
Theorem 5.3.1
Theorem 5.3.1. Let \(X_1, \ldots, X_n\) be a random sample from a \(N(\mu, \sigma^2)\) distribution. Then
- \(\bar{X}\) and \(S^2\) are independent random variables,
- \(\bar{X}\) has a \(N(\mu, \sigma^2/n)\) distribution,
- \((n-1)S^2/\sigma^2\) has a \(\chi^2_{n-1}\) distribution.
Part (b) was shown with the mgf in the previous section.
Reduction to \(\mu = 0\), \(\sigma^2 = 1\). Let \(Z_i = (X_i - \mu)/\sigma\). Then \(Z_1, \ldots, Z_n\) is a random sample from \(N(0,1)\), and \[ X_i = \mu + \sigma Z_i, \qquad \bar{X} = \mu + \sigma \bar{Z}, \qquad X_i - \bar{X} = \sigma(Z_i - \bar{Z}), \qquad S_X^2 = \sigma^2 S_Z^2. \] If \(\bar{Z}\) and \(S_Z^2\) are independent, so are \(\bar{X}\) and \(S_X^2\), since they are functions of \(\bar{Z}\) and \(S_Z^2\) respectively (Theorem 4.6.12). And \((n-1)S_X^2/\sigma^2 = (n-1)S_Z^2\). It is therefore enough to prove (a) and (c) for \(\mu = 0\) and \(\sigma^2 = 1\), which we assume from now on. CB refers to this as a location–scale argument (§3.5).
Proof of (a): independence of \(\bar{X}\) and \(S^2\)
Step 1: reduce. Since the deviations sum to zero, \(X_1 - \bar{X} = -\sum_{i=2}^{n}(X_i - \bar{X})\), and \[ S^2 = \frac{1}{n-1}\left(\left[\sum_{i=2}^{n}(X_i - \bar{X})\right]^2 + \sum_{i=2}^{n}(X_i - \bar{X})^2\right). \] So \(S^2\) is a function of \((X_2 - \bar{X}, \ldots, X_n - \bar{X})\) only. It is enough to show that this random vector is independent of \(\bar{X}\) (Theorem 4.6.12).
Step 2: transform. Define \[ y_1 = \bar{x}, \qquad y_k = x_k - \bar{x}, \quad k = 2, \ldots, n. \] The inverse transformation is \[ x_1 = y_1 - \sum_{i=2}^{n} y_i, \qquad x_k = y_k + y_1, \quad k = 2, \ldots, n, \] where the first equation uses \(x_1 - \bar{x} = -\sum_{i=2}^n (x_i - \bar{x})\). The Jacobian of the inverse transformation is \[ J = \begin{vmatrix} 1 & -1 & -1 & \cdots & -1 \\ 1 & 1 & 0 & \cdots & 0 \\ 1 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 0 & 0 & \cdots & 1 \end{vmatrix} = n. \] To see this, add rows \(2, \ldots, n\) to the first row. The first row becomes \((n, 0, \ldots, 0)\), and expanding along it leaves \(n\) times the determinant of the \((n-1) \times (n-1)\) identity matrix. CB states the Jacobian as \(1/n\): that is the determinant of the forward transformation from \(x\) to \(y\). The density formula uses the inverse transformation, whose determinant is \(n\).
Step 3: substitute. The joint pdf of the sample is \[ f(x_1, \ldots, x_n) = \frac{1}{(2\pi)^{n/2}} \exp\left(-\frac{1}{2}\sum_{i=1}^{n} x_i^2\right). \] In the new variables, \[ \begin{aligned} \sum_{i=1}^{n} x_i^2 &= \Bigl(y_1 - \sum_{i=2}^{n} y_i\Bigr)^2 + \sum_{i=2}^{n}(y_i + y_1)^2 \\ &= y_1^2 - 2y_1\sum_{i=2}^{n} y_i + \Bigl(\sum_{i=2}^{n} y_i\Bigr)^2 + \sum_{i=2}^{n} y_i^2 + 2y_1\sum_{i=2}^{n} y_i + (n-1)y_1^2 \\ &= n y_1^2 + \sum_{i=2}^{n} y_i^2 + \Bigl(\sum_{i=2}^{n} y_i\Bigr)^2. \end{aligned} \] Multiplying by the Jacobian \(n = n^{1/2} \cdot n^{1/2}\), \[ \begin{aligned} f(y_1, \ldots, y_n) = {} &\left[\left(\frac{n}{2\pi}\right)^{1/2} \exp\left(-\frac{n y_1^2}{2}\right)\right] \\ &\times \left[\frac{n^{1/2}}{(2\pi)^{(n-1)/2}} \exp\left(-\frac{1}{2}\left[\sum_{i=2}^{n} y_i^2 + \Big(\sum_{i=2}^{n} y_i\Big)^2\right]\right)\right]. \end{aligned} \]
Step 4: factor. The first bracket depends only on \(y_1 = \bar{x}\) and the second only on \((y_2, \ldots, y_n) = (x_2 - \bar{x}, \ldots, x_n - \bar{x})\). By Theorem 4.6.11, \(\bar{X}\) is independent of \((X_2 - \bar{X}, \ldots, X_n - \bar{X})\), and by Step 1, \(\bar{X}\) and \(S^2\) are independent. \(\qquad \blacksquare\)
The first bracket is the \(N(0, 1/n)\) density, which agrees with part (b).
Independence of \(\bar{X}\) and \(S^2\) is unexpected at first, since \(S^2\) is computed from \(\bar{X}\). It holds because of the form of the normal density. For other populations it fails: Geary’s theorem states that \(\bar{X}\) and \(S^2\) are independent if and only if the sample is iid normal.
Proof of (c): the distribution of \((n-1)S^2/\sigma^2\)
Write \(\bar{X}_k\) and \(S_k^2\) for the sample mean and sample variance of the first \(k\) observations. The ordering of the observations plays no role; any fixed ordering works. The proof uses the recursion \[ (n-1)S_n^2 = (n-2)S_{n-1}^2 + \frac{n-1}{n}\left(X_n - \bar{X}_{n-1}\right)^2. \tag{$\ast$} \]
Derivation of \((\ast)\). Let \(d = X_n - \bar{X}_{n-1}\). The mean satisfies \[ \bar{X}_n = \frac{(n-1)\bar{X}_{n-1} + X_n}{n} = \bar{X}_{n-1} + \frac{X_n - \bar{X}_{n-1}}{n} = \bar{X}_{n-1} + \frac{d}{n}. \] Split the sum of squares into the first \(n-1\) terms and the last one: \[ (n-1)S_n^2 = \sum_{i=1}^{n-1}(X_i - \bar{X}_n)^2 + (X_n - \bar{X}_n)^2. \] For the last term, \[ (X_n - \bar{X}_n)^2 = \left(X_n - \bar{X}_{n-1} - \frac{d}{n}\right)^2 = \left(d - \frac{d}{n}\right)^2 = \frac{(n-1)^2}{n^2}\, d^2. \] For the first \(n-1\) terms, \(X_i - \bar{X}_n = (X_i - \bar{X}_{n-1}) - d/n\), and the deviations \(X_i - \bar{X}_{n-1}\), \(i = 1, \ldots, n-1\), sum to zero: \[ \begin{aligned} \sum_{i=1}^{n-1}\left[(X_i - \bar{X}_{n-1}) - \frac{d}{n}\right]^2 &= \sum_{i=1}^{n-1}(X_i - \bar{X}_{n-1})^2 - \frac{2d}{n}\underbrace{\sum_{i=1}^{n-1}(X_i - \bar{X}_{n-1})}_{= 0} + (n-1)\frac{d^2}{n^2} \\ &= (n-2)S_{n-1}^2 + \frac{n-1}{n^2}\, d^2. \end{aligned} \] Adding the two parts, \[ (n-1)S_n^2 = (n-2)S_{n-1}^2 + \frac{(n-1) + (n-1)^2}{n^2}\, d^2 = (n-2)S_{n-1}^2 + \frac{n-1}{n}\, d^2, \] which is \((\ast)\).
Base case, \(n = 2\). With \(\bar{X}_2 = (X_1 + X_2)/2\), the two deviations are \(\pm (X_2 - X_1)/2\), so \[ S_2^2 = \frac{1}{2 - 1}\left(\frac{(X_2 - X_1)^2}{4} + \frac{(X_2 - X_1)^2}{4}\right) = \frac{1}{2}(X_2 - X_1)^2 = \left(\frac{X_2 - X_1}{\sqrt{2}}\right)^2. \] Since \(X_2 - X_1 \sim N(0, 2)\), the quantity \((X_2 - X_1)/\sqrt{2}\) is \(N(0,1)\), and \(S_2^2 \sim \chi^2_1\) by Lemma 5.3.2(a).
Inductive step. Assume that \((k-1)S_k^2 \sim \chi^2_{k-1}\). By \((\ast)\) with \(n = k+1\), \[ k S_{k+1}^2 = (k-1)S_k^2 + \frac{k}{k+1}\left(X_{k+1} - \bar{X}_k\right)^2. \]
The second term is \(\chi^2_1\). \(X_{k+1} \sim N(0,1)\) and \(\bar{X}_k \sim N(0, 1/k)\) are independent, so \(X_{k+1} - \bar{X}_k \sim N\bigl(0, 1 + \frac{1}{k}\bigr) = N\bigl(0, \frac{k+1}{k}\bigr)\). Standardizing and squaring, \(\frac{k}{k+1}(X_{k+1} - \bar{X}_k)^2 \sim \chi^2_1\).
The two terms are independent. By part (a) applied to the first \(k\) observations, \(S_k^2\) and \(\bar{X}_k\) are independent. \(X_{k+1}\) is independent of \(X_1, \ldots, X_k\), and so of \((\bar{X}_k, S_k^2)\). Hence \(S_k^2\), \(\bar{X}_k\) and \(X_{k+1}\) are mutually independent, and \(S_k^2\) is independent of \(X_{k+1} - \bar{X}_k\) (Theorem 4.6.12).
By Lemma 5.3.2(b), \(kS_{k+1}^2 \sim \chi^2_{k-1+1} = \chi^2_k\), which completes the induction. \(\qquad \blacksquare\)
The MSE of \(S^2\) and \(\hat\sigma^2\)
Part (c) gives the variance of \(S^2\). Since \(\operatorname{Var}(\chi^2_{n-1}) = 2(n-1)\), \[ \operatorname{Var}(S^2) = \frac{\sigma^4}{(n-1)^2}\operatorname{Var}\left(\frac{(n-1)S^2}{\sigma^2}\right) = \frac{\sigma^4}{(n-1)^2} \cdot 2(n-1) = \frac{2\sigma^4}{n-1}. \] \(S^2\) is unbiased, so \(\text{MSE}(S^2) = 2\sigma^4/(n-1)\).
Now consider the estimator with divisor \(n\), \[ \hat\sigma^2 = \frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X})^2 = \frac{n-1}{n}S^2. \] Its bias is \(E(\hat\sigma^2) - \sigma^2 = \frac{n-1}{n}\sigma^2 - \sigma^2 = -\sigma^2/n\), and its variance is \(\bigl(\frac{n-1}{n}\bigr)^2 \frac{2\sigma^4}{n-1} = \frac{2(n-1)}{n^2}\sigma^4\). By the bias–variance decomposition, \[ \text{MSE}(\hat\sigma^2) = \frac{2(n-1)}{n^2}\sigma^4 + \frac{\sigma^4}{n^2} = \frac{2n-1}{n^2}\,\sigma^4. \] Comparing, \[ \frac{2n-1}{n^2} < \frac{2}{n-1} \iff (2n-1)(n-1) < 2n^2 \iff 2n^2 - 3n + 1 < 2n^2 \iff 3n > 1, \] which holds for every \(n\). So \(\hat\sigma^2\) is biased, but has a smaller MSE than \(S^2\) for every sample size. \(\hat\sigma^2\) is also the maximum likelihood estimator of \(\sigma^2\) (Topic 6). Which of the two to use depends on the criterion; \(S^2\) is the standard choice because the exact distribution theory in this topic is stated in terms of it.
Why \(n - 1\)
The deviations from the sample mean always sum to zero, \[ \sum_{i=1}^{n}(X_i - \bar{X}) = 0, \] so once \(n-1\) of them are known, the last one is determined. Only \(n-1\) of the deviations are free, and the number \(n-1\) of degrees of freedom appears for this reason in three places.
In \(E(S^2) = \sigma^2\). For each \(i\), \(X_i - \bar{X} = \bigl(1 - \frac{1}{n}\bigr)X_i - \frac{1}{n}\sum_{j \neq i} X_j\), which has mean zero and variance \(\bigl(1 - \frac{1}{n}\bigr)^2\sigma^2 + \frac{n-1}{n^2}\sigma^2 = \frac{n-1}{n}\sigma^2\). Summing over \(i\) gives \(E\sum_{i=1}^{n}(X_i - \bar{X})^2 = (n-1)\sigma^2\). Measuring spread around \(\bar{X}\), which is fitted to the data, rather than around \(\mu\) makes the sum of squares smaller on average by one \(\sigma^2\).
In the proof of (a). \(S^2\) is a function of the \(n-1\) variables \(X_2 - \bar{X}, \ldots, X_n - \bar{X}\) alone, and these are independent of \(\bar{X}\).
In the proof of (c). The induction builds \((n-1)S^2/\sigma^2\) from \(n-1\) independent \(\chi^2_1\) terms, one for each observation after the first, so it is \(\chi^2_{n-1}\).
The \(t\)-distribution
If \(\sigma\) were known, \((\bar{X} - \mu)/(\sigma/\sqrt{n}) \sim N(0,1)\) could be used to judge how far \(\bar{X}\) is from \(\mu\). In practice \(\sigma\) is unknown and is replaced by \(S\). The resulting statistic is no longer normal.
Theorem (formulated as Definition 5.3.4 in CB). Let \(X_1, \ldots, X_n\) be a random sample from a \(N(\mu, \sigma^2)\) distribution. Then \[ T = \frac{\bar{X} - \mu}{S/\sqrt{n}} \] has a \(t\)-distribution with \(p = n-1\) degrees of freedom, with density \[ f_T(t) = \frac{\Gamma\!\left(\frac{p+1}{2}\right)} {\Gamma\!\left(\frac{p}{2}\right)\sqrt{p\pi}} \left(1 + \frac{t^2}{p}\right)^{-(p+1)/2}, \qquad -\infty < t < \infty. \]
Proof. Divide numerator and denominator by \(\sigma/\sqrt{n}\): \[ T = \frac{(\bar{X} - \mu)/(\sigma/\sqrt{n})}{\sqrt{S^2/\sigma^2}}. \] The numerator is \(N(0,1)\) by Theorem 5.3.1(b). The denominator is \(\sqrt{V/p}\) with \(V = (n-1)S^2/\sigma^2 \sim \chi^2_p\) by Theorem 5.3.1(c). Numerator and denominator are independent by Theorem 5.3.1(a). It remains to find the distribution of \(U/\sqrt{V/p}\), where \(U \sim N(0,1)\) and \(V \sim \chi^2_p\) are independent.
The joint pdf of \((U, V)\) is \[ f_{U,V}(u,v) = \frac{1}{(2\pi)^{1/2}} e^{-u^2/2}\, \frac{1}{\Gamma\!\left(\frac{p}{2}\right) 2^{p/2}}\, v^{p/2 - 1} e^{-v/2}, \qquad -\infty < u < \infty,\ 0 < v < \infty. \] Make the transformation \[ t = \frac{u}{\sqrt{v/p}}, \qquad w = v, \qquad\text{with inverse}\qquad u = t\left(\frac{w}{p}\right)^{1/2}, \qquad v = w. \] The Jacobian of the inverse is \[ \begin{vmatrix} \partial u/\partial t & \partial u/\partial w \\ \partial v/\partial t & \partial v/\partial w \end{vmatrix} = \begin{vmatrix} (w/p)^{1/2} & \ast \\ 0 & 1 \end{vmatrix} = \left(\frac{w}{p}\right)^{1/2}. \] The marginal pdf of \(T\) is obtained by integrating out \(w\). With \(u^2 = t^2 w/p\), \[ \begin{aligned} f_T(t) &= \int_0^{\infty} f_{U,V}\!\left(t\left(\tfrac{w}{p}\right)^{1/2}, w\right) \left(\tfrac{w}{p}\right)^{1/2} dw \\ &= \frac{1}{(2\pi)^{1/2}\,\Gamma\!\left(\frac{p}{2}\right) 2^{p/2}\, p^{1/2}} \int_0^{\infty} w^{\frac{p+1}{2} - 1} \exp\left(-\frac{1}{2}\left(1 + \frac{t^2}{p}\right)w\right) dw. \end{aligned} \] The integrand is a gamma kernel: \(\int_0^\infty w^{a-1}e^{-w/b}\,dw = \Gamma(a)\,b^a\), here with \(a = \frac{p+1}{2}\) and \(b = 2/(1 + t^2/p)\). Hence \[ f_T(t) = \frac{\Gamma\!\left(\frac{p+1}{2}\right)} {(2\pi)^{1/2}\,\Gamma\!\left(\frac{p}{2}\right) 2^{p/2}\, p^{1/2}} \left(\frac{2}{1 + t^2/p}\right)^{(p+1)/2}. \] The powers of 2 cancel, since \(2^{(p+1)/2}/(2^{1/2}\, 2^{p/2}) = 1\), and what is left is the stated density. \(\qquad \blacksquare\)
Some properties of the \(t_p\) distribution:
- It is symmetric around zero, with heavier tails than the normal.
- For \(p = 1\) the density is \(1/\bigl(\pi(1 + t^2)\bigr)\), the Cauchy distribution, which has no mean.
- \(E(T) = 0\) for \(p > 1\), and \(\operatorname{Var}(T) = p/(p-2)\) for \(p > 2\).
- As \(p \to \infty\), \((1 + t^2/p)^{-(p+1)/2} \to e^{-t^2/2}\), and \(t_p\) approaches \(N(0,1)\). For large samples, replacing \(\sigma\) by \(S\) makes little difference.
The \(F\)-distribution
Definition 5.3.6 (\(F\)-distribution). Let \(X_1, \ldots, X_n\) be a random sample from a \(N(\mu_X, \sigma_X^2)\) population, and let \(Y_1, \ldots, Y_m\) be a random sample from an independent \(N(\mu_Y, \sigma_Y^2)\) population. The random variable \[ F = \frac{S_X^2/\sigma_X^2}{S_Y^2/\sigma_Y^2} \] has an \(F\)-distribution with \((n-1, m-1)\) degrees of freedom. The pdf of an \(F_{p,q}\) variable is \[ f_F(x) = \frac{\Gamma\!\left(\frac{p+q}{2}\right)} {\Gamma\!\left(\frac{p}{2}\right)\Gamma\!\left(\frac{q}{2}\right)} \left(\frac{p}{q}\right)^{p/2} \frac{x^{p/2 - 1}}{\left(1 + \frac{p}{q}x\right)^{(p+q)/2}}, \qquad 0 < x < \infty. \]
By Theorem 5.3.1(c), \(S_X^2/\sigma_X^2\) is a \(\chi^2_{n-1}\) variable divided by \(n-1\), and \(S_Y^2/\sigma_Y^2\) is an independent \(\chi^2_{m-1}\) variable divided by \(m-1\). So \(F\) has the form \((U/p)/(V/q)\) with \(U \sim \chi^2_p\) and \(V \sim \chi^2_q\) independent, \(p = n-1\) and \(q = m-1\).
Derivation of the density. The argument is the same as for the \(t\)-density. Make the transformation \[ x = \frac{u/p}{v/q}, \qquad w = v, \qquad\text{with inverse}\qquad u = \frac{p}{q}\,x w, \qquad v = w, \] and Jacobian \((p/q)\,w\). Then \[ \begin{aligned} f_F(x) &= \int_0^{\infty} f_U\!\left(\tfrac{p}{q}xw\right) f_V(w)\, \tfrac{p}{q}\, w \; dw \\ &= \frac{(p/q)^{p/2}\, x^{p/2 - 1}} {\Gamma\!\left(\frac{p}{2}\right)\Gamma\!\left(\frac{q}{2}\right) 2^{(p+q)/2}} \int_0^{\infty} w^{\frac{p+q}{2} - 1} \exp\left(-\frac{1}{2}\left(1 + \frac{p}{q}x\right)w\right) dw. \end{aligned} \] The integral is again a gamma kernel, equal to \(\Gamma\!\left(\frac{p+q}{2}\right)\bigl(2/(1 + \frac{p}{q}x)\bigr)^{(p+q)/2}\). The powers of 2 cancel, which gives \(f_F\). \(\qquad \blacksquare\)
The representation \(F = (U/p)/(V/q)\) also gives two useful facts (CB, Theorem 5.3.8):
- If \(X \sim F_{p,q}\), then \(1/X \sim F_{q,p}\): the reciprocal swaps the roles of numerator and denominator.
- If \(T \sim t_q\), then \(T^2 \sim F_{1,q}\): write \(T = U/\sqrt{V/q}\) with \(U \sim N(0,1)\), so \(T^2 = (U^2/1)/(V/q)\) with \(U^2 \sim \chi^2_1\).
What these distributions are for
Each of the results above gives a function of the data and the parameter whose distribution is fully known.
Pivot (CB §9.2.2). A random variable \(Q = Q(X_1, \ldots, X_n, \theta)\) is a pivotal quantity, or pivot, if its distribution does not depend on any unknown parameter.
For a normal sample, \[ \frac{\bar{X} - \mu}{S/\sqrt{n}} \sim t_{n-1} \qquad\text{and}\qquad \frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1} \] are pivots for \(\mu\) and \(\sigma^2\): the right-hand distributions contain no unknown parameters. A probability statement about the pivot can be rearranged into a statement about the parameter.
Below, \(t_{p,\alpha}\) and \(\chi^2_{p,\alpha}\) denote upper quantiles: \(P(T > t_{p,\alpha}) = \alpha\) for \(T \sim t_p\), and similarly for \(\chi^2_p\).
Interval for \(\mu\) with \(\sigma\) unknown. Since \(T \sim t_{n-1}\) and the \(t\)-distribution is symmetric, \[ \begin{aligned} 1 - \alpha &= P\left(-t_{n-1,\alpha/2} \leq \frac{\bar{X} - \mu}{S/\sqrt{n}} \leq t_{n-1,\alpha/2}\right) \\ &= P\left(\bar{X} - t_{n-1,\alpha/2}\frac{S}{\sqrt{n}} \leq \mu \leq \bar{X} + t_{n-1,\alpha/2}\frac{S}{\sqrt{n}}\right). \end{aligned} \] The random interval \(\bar{X} \pm t_{n-1,\alpha/2}\, S/\sqrt{n}\) contains \(\mu\) with probability \(1 - \alpha\). This is the \(t\) interval.
Interval for \(\sigma^2\). Since \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\), \[ \begin{aligned} 1 - \alpha &= P\left(\chi^2_{n-1,1-\alpha/2} \leq \frac{(n-1)S^2}{\sigma^2} \leq \chi^2_{n-1,\alpha/2}\right) \\ &= P\left(\frac{(n-1)S^2}{\chi^2_{n-1,\alpha/2}} \leq \sigma^2 \leq \frac{(n-1)S^2}{\chi^2_{n-1,1-\alpha/2}}\right). \end{aligned} \] Taking reciprocals reverses the inequalities, so the larger quantile gives the lower endpoint. The \(\chi^2\) distribution is not symmetric, so the interval is not symmetric around \(S^2\).
Example. A sample of \(n = 10\) observations from a normal population gives \(\bar{x} = 5.2\) and \(s = 1.5\). From tables, \(t_{9, 0.025} = 2.262\), \(\chi^2_{9, 0.025} = 19.023\) and \(\chi^2_{9, 0.975} = 2.700\).
A 95% interval for \(\mu\) is \[ 5.2 \pm 2.262 \cdot \frac{1.5}{\sqrt{10}} = 5.2 \pm 1.07 = [4.13,\ 6.27]. \]
With \((n-1)s^2 = 9 \cdot 2.25 = 20.25\), a 95% interval for \(\sigma^2\) is \[ \left[\frac{20.25}{19.023},\ \frac{20.25}{2.700}\right] = [1.06,\ 7.50], \] and taking square roots, \([1.03,\ 2.74]\) for \(\sigma\).
Test of \(\sigma_X^2 = \sigma_Y^2\). For two independent normal samples, the ratio \(F = (S_X^2/\sigma_X^2)/(S_Y^2/\sigma_Y^2)\) contains the unknown variances. Under the hypothesis \(\sigma_X^2 = \sigma_Y^2\) they cancel, and \(S_X^2/S_Y^2 \sim F_{n-1,m-1}\) can be computed from the data alone. A value far out in either tail of \(F_{n-1,m-1}\) is evidence against the hypothesis. Hypothesis testing is treated in Topic 8.
Order statistics
Every statistic so far has been an average, or built from averages. The sample mean estimates the population mean, the sample variance estimates the population variance, and the \(t\)-statistic compares the two. Some questions need a different kind of statistic.
The German tank problem. During the Second World War, the Allies wanted to estimate the number of tanks Germany produced. Captured and destroyed tanks carried serial numbers. Suppose the numbers are \(1, 2, \ldots, N\), that the observed tanks are a random sample drawn without replacement, and that the observed serial numbers are \(19, 40, 42\) and \(60\). How should \(N\) be estimated?
An average of the four numbers is not a sensible basis. Whatever the first three values are, \(N\) is at least 60. The natural statistic is the maximum.
Definition 5.4.1. The order statistics of a random sample \(X_1, \ldots, X_n\) are the sample values placed in ascending order, denoted \(X_{(1)} \leq X_{(2)} \leq \cdots \leq X_{(n)}\). In particular \(X_{(1)} = \min_i X_i\) and \(X_{(n)} = \max_i X_i\).
The German tank estimator
Draws without replacement are not independent, so the tank sample is not iid. The distribution of the maximum can still be found by counting. All \(\binom{N}{n}\) subsets of size \(n\) are equally likely. The maximum equals \(x\) exactly when the subset contains \(x\) and \(n-1\) numbers from \(\{1, \ldots, x-1\}\), so \[ P(X_{(n)} = x) = \frac{\binom{x-1}{n-1}}{\binom{N}{n}}, \qquad x = n, \ldots, N. \]
For the expectation, use \(x\binom{x-1}{n-1} = n\binom{x}{n}\) and the identity \(\sum_{x=n}^{N}\binom{x}{n} = \binom{N+1}{n+1}\): \[ E X_{(n)} = \frac{1}{\binom{N}{n}}\sum_{x=n}^{N} x\binom{x-1}{n-1} = \frac{n}{\binom{N}{n}}\sum_{x=n}^{N}\binom{x}{n} = \frac{n\binom{N+1}{n+1}}{\binom{N}{n}} = \frac{n(N+1)}{n+1}. \] The identity counts the subsets of size \(n+1\) of \(\{1, \ldots, N+1\}\) according to their largest element \(x + 1\).
Solving \(E X_{(n)} = n(N+1)/(n+1)\) for \(N\) suggests the estimator \[ \hat{N} = \frac{n+1}{n}\, X_{(n)} - 1, \] which is unbiased by construction: \(E\hat{N} = \frac{n+1}{n}\cdot\frac{n(N+1)}{n+1} - 1 = N\). With \(n = 4\) and observed maximum \(60\), \[ \hat{N} = \frac{5}{4} \cdot 60 - 1 = 74. \]
The estimator can be read as the observed maximum plus the average gap between observed serial numbers: \(\hat{N} = m + (m - n)/n\), where \(m - n\) is the number of unobserved serial numbers below \(m\).
For June 1940 to September 1942, conventional intelligence estimated German tank production at about 1,400 a month. The serial-number estimate was 246. German records captured after the war showed 245.
Continuous populations
For a random sample (iid, so with replacement) from a continuous population with cdf \(F\) and pdf \(f\), the distribution of the maximum follows from independence: \[ F_{X_{(n)}}(x) = P(X_{(n)} \leq x) = P(X_1 \leq x, \ldots, X_n \leq x) = \bigl(F(x)\bigr)^n, \] and differentiating, \[ f_{X_{(n)}}(x) = n f(x) \bigl(F(x)\bigr)^{n-1}. \] The minimum is treated in the same way in problem P10. The general result is:
Theorem 5.4.4. Let \(X_{(1)}, \ldots, X_{(n)}\) be the order statistics of a random sample from a continuous population with cdf \(F\) and pdf \(f\). The pdf of \(X_{(j)}\) is \[ f_{X_{(j)}}(x) = \frac{n!}{(j-1)!\,(n-j)!}\, f(x)\, \bigl(F(x)\bigr)^{j-1} \bigl(1 - F(x)\bigr)^{n-j}. \]
For \(j = n\) this is the density of the maximum above. The argument behind the general formula is problem A3.