Probability · from the model to the data
Statistics · from the data to the model
So far, we have been analyzing probability models and asking: given a distribution, what are the probabilities of various outcomes?
\[\longrightarrow \text{Let us call this \textbf{probability theory}.}\]
We will now reverse the direction. We observe data, and from that we will infer properties of the unknown probability model that generated the data.
\[\longrightarrow \text{Let us call this \textbf{statistics}, or \textbf{statistical inference}.}\]
We will concentrate on parametric statistics: we are willing to assume that the observations come from a specific probability model that is known up to a set of parameters. The task is then to estimate the parameters using data.
Examples:
Other times we are not willing to make such assumptions and try to estimate \(f(\cdot)\) directly from data — that is called nonparametric statistics.
Definition 7.1.1. A point estimator is any function \(W(X_1, \ldots, X_n)\) of a sample. That is, any statistic is a point estimator.
This definition is very general — it does not reference what \(W(\cdot)\) is an estimator of. It is simply a function of the sample.
Two important distinctions:
The most widely used strategy is maximum likelihood. Two simpler strategies serve as useful illustrations.
1. Use logic. We already know that the sample mean converges to the population mean, both strongly and weakly. If \(X_1, \ldots, X_n \overset{\text{iid}}{\sim} N(\mu, \sigma^2)\), then \(E(X_i) = \mu\), and hence \(\bar{X} \xrightarrow{\text{a.s.}} \mu\). We also know that \(E(\bar{X}) = \mu\). The sample mean is therefore a very natural estimator of \(\mu\).
2. The method of moments. Equate the population mean and the sample mean — the theory we have learnt tells us these should be close. The same body of theory provides the same results for higher moments: \[m_k = \frac{1}{n}\sum_{i=1}^n X_i^k \approx \mu_k = E X^k.\]
If we have \(p\) unknown parameters, we write up the system of equations \[m_k = \mu_k(\theta_1, \ldots, \theta_p), \quad k = 1, \ldots, p,\] and solve for \(\theta_k\), \(k = 1, \ldots, p\).
Example 7.2.1. Suppose \(X_1, \ldots, X_n\) are iid \(N(\theta, \sigma^2)\). The two empirical moments are \(m_1 = \bar{X}\) and \(m_2 = \tfrac{1}{n}\sum X_i^2\). The population moments are \(\mu_1 = E(X) = \theta\) and \(\mu_2 = E(X^2) = \theta^2 + \sigma^2\).
We must therefore solve the system \(\bar{X} = \theta\) and \(\tfrac{1}{n}\sum X_i^2 = \theta^2 + \sigma^2\), yielding the method of moments estimators: \[\tilde{\theta} = \bar{X} \qquad \text{and} \qquad \tilde{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n X_i^2 - \bar{X}^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2.\]
Suppose we toss a coin 10 times and observe heads on 8 of them. We model the tosses as iid Bernoulli\((\theta)\) and ask: which value of \(\theta\) best explains what we just saw?
For any \(\theta\), the probability of the observed event is \[P_\theta(X = 8) = \binom{10}{8}\theta^8(1-\theta)^2.\]
| \(\theta\) | \(P_\theta(X = 8)\) |
|---|---|
| 0.5 | 0.044 |
| 0.7 | 0.233 |
| 0.8 | 0.302 |
| 0.9 | 0.194 |
| 1.0 | 0 |
The data are fixed and the parameter is varying. Among these values, \(\theta = 0.8\) makes the data most probable. It can be shown that the maximum of \(\theta \mapsto \theta^8(1-\theta)^2\) is achieved exactly at \(\theta = 0.8\).
The reasoning above is general. For any fixed dataset, the map \[\theta \mapsto P_\theta(\text{observed data})\] gives an objective ranking of candidate values of \(\theta\): those that make what we actually saw more probable are preferred over those that make it less probable. The model \(\{P_\theta\}_{\theta \in \Theta}\) provides the ranking — we do not have to introduce any prior preference of our own.
Definition (Likelihood, discrete case). Let \(X\) have probability mass function \(p_\theta(x)\), indexed by \(\theta \in \Theta\). Given an observed value \(x\), the likelihood function is \[L(\theta;\, x) = p_\theta(x), \quad \theta \in \Theta,\] viewed as a function of \(\theta\) for the fixed data \(x\).
We write just \(L(\theta)\) when the data are understood. The likelihood is not a probability distribution over \(\Theta\) and is not normalised. But the ratio \(L(\theta_2)/L(\theta_1)\) has a clear meaning: it tells us how many times more probable the data are under \(\theta_2\) than under \(\theta_1\).
If \(X\) has density \(p_\theta(x)\), then \(P_\theta(X = x) = 0\) and the definition above appears to collapse. The fix is to replace the point \(x\) by a small interval: \[P(X \in [x-h, x+h]) \approx 2h\, p_\theta(x).\]
The factor \(2h\) does not depend on \(\theta\), so for comparing values of \(\theta\) it is irrelevant — it cancels in every likelihood ratio. We therefore define the continuous likelihood by the same formula \(L(\theta;\, x) = p_\theta(x)\), understanding likelihood only up to a constant that does not involve \(\theta\).
For an iid sample \(X_1, \ldots, X_n\) with a common density (or mass) \(p_\theta(x)\), independence makes the likelihood a product: \[L(\theta;\, X_1, \ldots, X_n) = \prod_{i=1}^n p_\theta(X_i).\]
We almost always work with its logarithm, which is a monotone function that does not change our order of preference, leading to the log-likelihood: \[\ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log p_\theta(X_i).\]
Products of small numbers turn into sums of manageable ones, and the quantities that govern inference — scores and information — are derivatives of \(\ell\), not \(L\).
Once we have \(\ell(\theta)\), the data are essentially summarised. Everything we will derive — the point estimate (where is \(\ell\) maximised?), the standard error (how sharply is \(\ell\) peaked there?), and likelihood ratios — can be read off the shape of \(\ell\).
The log-likelihood is a package containing all the information in the data about the unknown parameter, given the model.
Definition (Maximum Likelihood Estimator). The maximum likelihood estimator (MLE) of \(\theta\) is any value \(\hat\theta \in \Theta\) such that \[\hat\theta \in \operatorname{arg\,max}_{\theta \in \Theta}\; \ell(\theta).\]
In the interior of the parameter space and under smoothness conditions, \(\hat\theta\) solves the score equation \[\ell'(\hat\theta) = 0,\] with the second-order check that the curvature is negative there.
The MLE is the natural endpoint of the ranking we set up at the start: among all candidate \(\theta\), we pick the one that makes the data most probable.
But the MLE is not the only thing the likelihood tells us. Expand \(\ell\) around \(\hat\theta\) to second order: \[\ell(\theta) \approx \ell(\hat\theta) - \tfrac{1}{2}\bigl(-\ell''(\hat\theta)\bigr)(\theta - \hat\theta)^2.\]
The first-order term vanishes because \(\ell'(\hat\theta) = 0\), and the remaining quantity \(-\ell''(\hat\theta)\) controls how quickly \(\ell\) falls away from its maximum:
The curvature at \(\hat\theta\) is therefore exactly the right object to measure how informative the data are about \(\theta\).
Let \(X_1, \ldots, X_n\) be iid \(N(\mu, \sigma^2)\) with \(\sigma^2\) known. The log-likelihood is \[\begin{align} \ell(\mu) &= \sum_{i=1}^n \log\!\left[\frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\left\{-\tfrac{1}{2}\left(\frac{X_i - \mu}{\sigma}\right)^2\right\}\right] \\ &= -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (X_i - \mu)^2. \end{align}\]
Differentiating: \(\ell'(\mu) = \dfrac{1}{\sigma^2}\displaystyle\sum_{i=1}^n (X_i - \mu)\).
Setting \(\ell'(\mu) = 0\) gives the MLE \(\hat\mu = \bar{X}_n\), the sample mean.
Definition (Score and observed information). The score function is \[U(\theta) = \ell'(\theta) = \frac{\partial}{\partial\theta}\log L(\theta).\]
The observed Fisher information is \[I(\theta) = -\ell''(\theta) = -\frac{\partial^2}{\partial\theta^2}\log L(\theta).\]
The MLE solves \(U(\hat\theta) = 0\). The number \(I(\hat\theta)\) is the curvature of \(\ell\) at its maximum.
There is also a model-level companion: the expected (or plain Fisher) information \[\mathcal{I}(\theta) = E_\theta\!\bigl(I(\theta)\bigr) \equiv \operatorname{Var}(U(\theta)),\] which averages the curvature over what we might have observed.
For inference, we use the observed information at the MLE, because it reflects the curvature of the log-likelihood that we actually computed.
We have alluded many times to the role that information plays when measuring precision of estimates. The link follows from two facts.
Fact 1: Variance of score = expected information (stated for a single observation, \(n = 1\)).
Starting from \(\int p_\theta(x)\,dx = 1\), differentiate under the integral: \[0 = \frac{d}{d\theta}\int p_\theta\,dx = \int p_\theta \frac{\partial \log p_\theta}{\partial\theta}\,dx = E_\theta(U(\theta)). \quad (*)\]
Differentiate once more: \[0 = \frac{d}{d\theta}\int p_\theta U\,dx = \int \frac{dp_\theta}{d\theta} U\,dx + \int p_\theta \frac{dU}{d\theta}\,dx = E_\theta(U^2) - E_\theta(I(\theta)),\]
so \(E_\theta(U^2) = E_\theta(I(\theta))\), and because \(E_\theta(U) = 0\) the left side is just \(\operatorname{Var}_\theta(U)\): \[\operatorname{Var}_\theta(U(\theta)) = E_\theta(I(\theta)) = \mathcal{I}(\theta).\]
(This argument assumes we can interchange the order of integration and differentiation — see Theorem 2.4.3 in CB for details, but this holds in most practical situations.)
Fact 2: The score satisfies a central limit theorem.
Let \(\theta_0\) denote the true parameter value. When we evaluate the score at \(\theta_0\), \[U(\theta_0) = \sum_{i=1}^n \frac{\partial}{\partial\theta}\log p_\theta(X_i)\big|_{\theta=\theta_0},\] it is a sum of \(n\) iid terms. By Fact 1 applied with \(\theta = \theta_0\), each term has mean zero and variance \(\mathcal{I}(\theta_0)\), so by the CLT: \[U(\theta_0)/\sqrt{n} \xrightarrow{d} N(0,\, \mathcal{I}(\theta_0)).\]
The MLE solves \(U(\hat\theta) = 0\). Taylor-expanding \(U\) around \(\theta_0\): \[0 = U(\hat\theta) \approx U(\theta_0) - I(\theta_0)(\hat\theta - \theta_0) \implies \hat\theta - \theta_0 \approx \frac{U(\theta_0)}{I(\theta_0)}.\]
The numerator has variance \(\approx n\mathcal{I}(\theta_0)\) (Fact 2); the denominator converges to \(n\mathcal{I}(\theta_0)\) by the LLN. Substituting: \[\operatorname{Var}(\hat\theta) \approx \frac{n\mathcal{I}(\theta_0)}{(n\mathcal{I}(\theta_0))^2} = \frac{1}{n\mathcal{I}(\theta_0)} \approx I(\hat\theta)^{-1}, \qquad \text{so}\quad \operatorname{SE}(\hat\theta) \approx \frac{1}{\sqrt{I(\hat\theta)}}.\]
Continuing the example above: \[\ell''(\mu) = -\frac{n}{\sigma^2}, \qquad I(\hat\mu) = -\ell''(\hat\mu) = \frac{n}{\sigma^2}.\]
Therefore: \[\operatorname{SE}(\hat\mu) \approx \frac{1}{\sqrt{I(\hat\mu)}} = \frac{\sigma}{\sqrt{n}},\]
which is exactly the standard error we already know for the sample mean.
Note the convergence rate: \(\hat\mu\) converges to \(\mu\) at speed \(1/\sqrt{n}\) — the same rate as the CLT predicts.
The whole construction generalizes naturally to a parameter vector \(\theta = (\theta_1, \ldots, \theta_p)\).
The score vector is the gradient: \[U(\theta) = \nabla\ell(\theta) = \left(\frac{\partial\ell}{\partial\theta_1}, \ldots, \frac{\partial\ell}{\partial\theta_p}\right) \in \mathbb{R}^p.\]
The observed information matrix is the negative Hessian matrix: \[I(\theta) = -\nabla^2\ell(\theta) = \left\{-\frac{\partial^2\ell}{\partial\theta_j\,\partial\theta_k}\right\}_{1 \leq j,k \leq p} \in \mathbb{R}^{p \times p}.\]
The MLE solves \(U(\hat\theta) = \mathbf{0}\), and the multivariate analogue of inverting the curvature gives: \[\widehat{\operatorname{Cov}}(\hat\theta) \approx I(\hat\theta)^{-1}.\]
The diagonal of \(I(\hat\theta)^{-1}\) gives the variances of individual parameter estimates \(\hat\theta_j\), and the off-diagonal elements provide the approximate covariances between parameter estimates.
Let \(X_1, \ldots, X_n\) be iid \(N(\mu, \sigma^2)\) with both \(\mu\) and \(\sigma^2\) unknown. Write \(\tau = \sigma^2\) to simplify notation when differentiating. The log-likelihood is \[\ell(\mu, \tau) = -\frac{n}{2}\log(2\pi) - \frac{n}{2}\log\tau - \frac{1}{2\tau}\sum_{i=1}^n(X_i - \mu)^2.\]
Differentiating and setting partial derivatives to zero: \[\begin{align} \frac{\partial\ell}{\partial\mu} &= \frac{1}{\tau}\sum_{i=1}^n(X_i - \mu) = 0 &&\Rightarrow \hat\mu = \bar{X}_n, \\ \frac{\partial\ell}{\partial\tau} &= -\frac{n}{2\tau} + \frac{1}{2\tau^2}\sum_{i=1}^n(X_i - \mu)^2 = 0 &&\Rightarrow \hat\sigma^2 = \frac{1}{n}\sum_{i=1}^n(X_i - \bar{X}_n)^2. \end{align}\]
Note the divisor is \(n\), not \(n-1\): maximum likelihood estimates are not always unbiased. They are, however, asymptotically unbiased under regularity conditions (the difference between \(1/n\) and \(1/(n-1)\) goes to zero as \(n \to \infty\)).
The second partial derivatives are: \[\frac{\partial^2\ell}{\partial\mu^2} = -\frac{n}{\tau}, \qquad \frac{\partial^2\ell}{\partial\mu\,\partial\tau} = -\frac{1}{\tau^2}\sum(X_i - \mu), \qquad \frac{\partial^2\ell}{\partial\tau^2} = \frac{n}{2\tau^2} - \frac{1}{\tau^3}\sum(X_i - \mu)^2.\]
At the MLE, \(\sum(X_i - \hat\mu) = 0\), making the off-diagonal element vanish. Using \(\sum(X_i - \bar{X}_n)^2 = n\hat\sigma^2\): \[I(\hat\theta) = \begin{pmatrix} n/\hat\sigma^2 & 0 \\ 0 & n/(2\hat\sigma^4) \end{pmatrix}, \qquad I(\hat\theta)^{-1} = \begin{pmatrix} \hat\sigma^2/n & 0 \\ 0 & 2\hat\sigma^4/n \end{pmatrix}.\]
Reading off the diagonal: \[\operatorname{SE}(\hat\mu) \approx \frac{\hat\sigma}{\sqrt{n}}, \qquad \operatorname{SE}(\hat\sigma^2) \approx \hat\sigma^2\sqrt{\frac{2}{n}}.\]
The off-diagonal zero is not surprising, given the classical result that the sample mean and sample variance are statistically independent.
Apply the likelihood machinery to ordinary simple linear regression, with residuals drawn from the normal distribution. We observe pairs \((x_i, y_i)\), \(i = 1, \ldots, n\), modelled as \[Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i, \qquad \varepsilon_i \overset{\text{iid}}{\sim} N(0, \sigma^2),\] with the \(x_i\) regarded as fixed. The unknown parameters are \(\theta = (\beta_0, \beta_1, \sigma^2)\), and the log-likelihood is \[\ell(\beta_0, \beta_1, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n(Y_i - \beta_0 - \beta_1 X_i)^2.\]
For a fixed \(\sigma^2\), maximizing over \((\beta_0, \beta_1)\) is the same as minimizing the sum of squared residuals: \[\operatorname{SSR}(\beta_0, \beta_1) = \sum_{i=1}^n(Y_i - \beta_0 - \beta_1 X_i)^2.\]
The MLE of the regression coefficients is the ordinary least squares estimator.
Setting \(\partial\ell/\partial\beta_0 = 0\) and \(\partial\ell/\partial\beta_1 = 0\) gives the normal equations: \[\sum_{i=1}^n(Y_i - \beta_0 - \beta_1 X_i) = 0 \qquad \text{and} \qquad \sum_{i=1}^n X_i(Y_i - \beta_0 - \beta_1 X_i) = 0.\]
Writing \(\bar{x} = n^{-1}\sum x_i\), \(\bar{y} = n^{-1}\sum y_i\), \(S_{xx} = \sum(x_i - \bar{x})^2\), \(S_{xy} = \sum(x_i - \bar{x})(y_i - \bar{y})\), solving the \(2 \times 2\) system gives the classical closed-form estimators: \[\hat\beta_1 = \frac{S_{xy}}{S_{xx}}, \qquad \hat\beta_0 = \bar{y} - \hat\beta_1\bar{x}.\]
Finally, plugging \((\hat\beta_0, \hat\beta_1)\) into \(\partial\ell/\partial\sigma^2 = 0\) gives: \[\hat\sigma^2 = \frac{1}{n}\sum_{i=1}^n(Y_i - \hat\beta_0 - \hat\beta_1 X_i)^2 = \frac{\operatorname{SSR}(\hat\beta_0, \hat\beta_1)}{n}.\]
A pleasant surprise: the MLEs of the regression coefficients under normality are the same as the least squares estimates, which can be derived without assuming anything about the distribution of the residuals!
Theorem (Asymptotic normality of the MLE). Let \(X_1, \ldots, X_n\) be iid with density \(f_\theta\), where \(\theta_0\) lies in the interior of \(\Theta \subset \mathbb{R}\). Assume the following regularity conditions: the support of \(f_\theta\) does not depend on \(\theta\); \(\log f_\theta\) is three times continuously differentiable in a neighbourhood of \(\theta_0\) with bounds permitting two interchanges of differentiation and integration; the per-observation Fisher information \(\mathcal{I}(\theta_0)\) is finite and positive; and the MLE is consistent: \(\hat\theta_n \xrightarrow{p} \theta_0\). Then: \[\sqrt{n}\,(\hat\theta_n - \theta_0) \xrightarrow{d} N\!\bigl(0,\, \mathcal{I}(\theta_0)^{-1}\bigr).\]
Apply Taylor’s theorem to the score \(U\) around \(\theta_0\). There exists \(\theta_n^*\) between \(\hat\theta_n\) and \(\theta_0\) such that \[0 = U(\hat\theta_n) = U(\theta_0) + U'(\theta_n^*)(\hat\theta_n - \theta_0).\]
Since \(U' = \ell'' = -I\), rearranging gives \(\hat\theta_n - \theta_0 = U(\theta_0)/I(\theta_n^*)\), and multiplying by \(\sqrt{n}\): \[\sqrt{n}\,(\hat\theta_n - \theta_0) = \frac{U(\theta_0)/\sqrt{n}}{I(\theta_n^*)/n}.\]
Numerator: \(U(\theta_0) = \sum_{i=1}^n \frac{\partial}{\partial\theta}\log p_\theta(X_i)\big|_{\theta=\theta_0}\) is a sum of \(n\) iid terms with mean zero and variance \(\mathcal{I}(\theta_0)\). By the CLT: \(U(\theta_0)/\sqrt{n} \xrightarrow{d} N(0,\, \mathcal{I}(\theta_0))\).
Denominator: \(I(\theta_n^*)/n\) is a sample average evaluated at \(\theta_n^*\), which is squeezed between \(\hat\theta_n\) and \(\theta_0\). Consistency of \(\hat\theta_n\) and continuity of \(\partial^2\log f_\theta/\partial\theta^2\) together with the LLN give \(I(\theta_n^*)/n \xrightarrow{p} \mathcal{I}(\theta_0)\).
Combine with Slutsky: \(\sqrt{n}(\hat\theta_n - \theta_0) \xrightarrow{d} N(0, \mathcal{I}(\theta_0))/\mathcal{I}(\theta_0) = N(0, \mathcal{I}(\theta_0)^{-1})\). \(\blacksquare\)
The proof relied on \(\hat\theta_n \xrightarrow{p} \theta_0\). Consistency is not automatic, but follows from a nice identity. By the LLN, the rescaled sample log-likelihood converges pointwise: \[\frac{1}{n}\ell_n(\theta) \xrightarrow{p} M(\theta) = E_{\theta_0}\!\bigl[\log f_\theta(X)\bigr].\]
The difference between \(M\) at the true value \(\theta_0\) and any other \(\theta\) is \[M(\theta_0) - M(\theta) = E_{\theta_0}\!\left[\log\frac{f_{\theta_0}(X)}{f_\theta(X)}\right] = D_{\text{KL}}(f_{\theta_0},\, f_\theta) \geq 0,\] with equality if and only if \(\theta = \theta_0\) almost surely. The non-negativity is Jensen’s inequality applied to \(-\log(\cdot)\).
Under identifiability (distinct \(\theta\) giving distinct distributions), \(M(\theta)\) is therefore strictly maximised at the true value \(\theta_0\). The sample log-likelihood approximates \(M\), and the sample maximizer \(\hat\theta_n\) must track the population maximizer \(\theta_0\).
The asymptotic normality theorem can be restated as the approximation \[\hat\theta_n \;\dot\sim\; N\!\left(\theta_0,\; \frac{1}{n\mathcal{I}(\theta_0)}\right).\]
Replacing \(n\mathcal{I}(\theta_0)\) by the observed information \(I(\hat\theta_n)\) (by the LLN and consistency/continuous mapping) and \(\theta_0\) by \(\hat\theta_n\), gives the working approximation: \[\hat\theta_n \;\dot\sim\; N\!\bigl(\theta_0,\; I(\hat\theta_n)^{-1}\bigr).\]
This is the same statement we used heuristically in the examples above, now justified as a limit theorem about the sampling distribution of the MLE. The classical confidence interval \[\hat\theta_n \pm z_{1-\alpha/2}\,\frac{1}{\sqrt{I(\hat\theta_n)}}\] falls out immediately and is the workhorse of large-sample inference.
Multivariate version. With \(\theta_0 \in \mathbb{R}^p\) in the interior of \(\Theta\), an entirely analogous argument using the multivariate CLT for the score vector \(U(\theta_0)\) and convergence of the information matrix \(I(\theta_n^*)/n\) to the per-observation information matrix \(\mathcal{I}(\theta_0)\) gives: \[\sqrt{n}\,(\hat\theta_n - \theta_0) \xrightarrow{d} N_p\!\bigl(\mathbf{0},\; \mathcal{I}(\theta_0)^{-1}\bigr).\]
The two-parameter normal example and the regression example above are both special cases.
Extension via the delta method. A smooth transformation of the MLE inherits asymptotic normality automatically. If \(g : \mathbb{R} \to \mathbb{R}\) is differentiable at \(\theta_0\) with \(g'(\theta_0) \neq 0\), the delta method gives: \[\sqrt{n}\,\bigl(g(\hat\theta_n) - g(\theta_0)\bigr) \xrightarrow{d} N\!\left(0,\; \frac{g'(\theta_0)^2}{\mathcal{I}(\theta_0)}\right).\]
We immediately obtain the limit distribution of any smooth transformation of the MLE: odds ratios, log-odds, predicted probabilities, or anything else we might want to report.
Theorem 7.2.10 (CB). If \(\hat\theta\) is the MLE of \(\theta\), then \(\tau(\hat\theta)\) is the MLE of \(\tau(\theta)\).
This result is not very surprising given the delta method above, but it is useful: it means we never have to re-derive the MLE for a transformed parameter. We can simply plug the MLE into the transformation.
Combined with the delta method, this gives us both the point estimate and the standard error of any smooth function of the parameters — directly from the MLE and the observed information matrix we have already computed.