6 · Maximum Likelihood Estimation
So far, we have been analyzing probability models and how we can say something about the probability of various outcomes if they come from this or that distribution.
\(\rightarrow\) Let us call this 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.
\(\rightarrow\) Let us call this statistics or statistical inference.
One more important remark: we will concentrate on parametric statistics, where 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. For example:
- Assume \(X_1, X_2, \ldots, X_n \sim\) iid \(N(\mu, \sigma^2)\). Estimate \(\mu\) and \(\sigma\).
- Assume \(Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i\), where \(\varepsilon_i \sim N(0, \sigma^2)\). Estimate \((\beta_0, \beta_1, \sigma)\) based on a sample \((Y_i, X_i)\), \(i = 1, \ldots, n\).
Other times we are not willing to make such assumptions, and we try to estimate \(f(\cdot)\) directly from data, e.g. if \(X_i \sim\) iid \(f(x)\) or if \(Y_i = f(X_i, \varepsilon_i)\), \(i = 1, \ldots, n\). 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, especially because it does not reference what \(W(\cdot)\) is an estimator of. It is simply a function of the sample.
An estimator is a function of a sample, and is therefore a random variable itself, having a distribution and other stochastic properties such as mean and variance.
An estimate is a realized value \(W(x_1, \ldots, x_n)\) of the estimator.
The most widely used strategy for finding estimators of parameters is to use maximum likelihood. Let us (as CB) mention two other strategies as 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 \sim\) iid \(N(\mu, \sigma^2)\), then we know that \(E(X_i) = \mu\), and hence that \(\bar{X} \overset{a.s.}{\to} \mu\). We also know that \(E(\bar{X}) = \mu\). The sample mean is therefore a very natural estimator of the population mean.
2. The method of moments. Continue the reasoning from the previous point where we equate the population mean and the sample mean because the theory we have learnt so far tells us that these numbers should be close to each other. The same body of theory provides the same results for higher moments; that is \[ m_k = \frac{1}{n}\sum_{i=1}^{n} X_i^k \approx \mu_k = E X^k. \]
The left hand side can be calculated from the data. The right hand side will be expressions containing the unknown parameters. If we have \(p\) unknown parameters, we can write up the system of equations \[ m_k = \mu_k(\theta_1, \ldots, \theta_p), \qquad 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)\). There are two unknown parameters in this model, and the two first empirical moments are \(m_1 = \bar{X}\) and \(m_2 = (1/n)\sum_{i=1}^n X_i^2\). The two first population moments are \(\mu_1 = E(X) = \theta\) and \(\mu_2 = E(X^2) = \theta^2 + \sigma^2\) (because \(\sigma^2 = E(X^2) - (E(X))^2\)).
We must therefore solve the system \[ \bar{X} = \theta, \qquad \frac{1}{n}\sum_{i=1}^{n} X_i^2 = \theta^2 + \sigma^2, \] which yields 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. \]
Likelihood — a motivating example
(The following is loosely based on Pawitan’s “In All Likelihood”.)
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. \]
Plugging in a few candidate values:
\[ \begin{array}{c|c} \theta & P_\theta(X = 8) \\ \hline 0.5 & 0.044 \\ 0.7 & 0.233 \\ 0.8 & 0.302 \\ 0.9 & 0.194 \\ 1.0 & 0 \end{array} \]
The data are fixed and the parameter is varying, with the probability of observing what we have observed varying with it. Among the values in the table, \(\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 that we used above is general. For any fixed dataset, the map \[ \theta \longrightarrow P_\theta(\text{Observed data}) \] gives us an objective ranking of the 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.
This map has a name:
The likelihood function
Definition (likelihood, the 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), \qquad \theta \in \Theta, \] viewed as a function of \(\theta\) for the fixed data \(x\).
We will write just \(L(\theta)\) when the data are understood. The likelihood is not a probability distribution over \(\theta\), and it 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\).
A remark about continuous models. 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\big(X \in [x-h,\, x+h]\big) \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\).
Log-likelihood of a sample. 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 (see below) — are derivatives of \(\ell\), not \(L\).
Information content of \(\ell(\theta)\). Once we have \(\ell(\theta)\), the data are essentially summarised. Everything we will derive below, such as 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.
Maximum likelihood estimation
Definition (maximum likelihood estimator). The maximum likelihood estimator (MLE) of \(\theta\) is any value \(\hat{\theta} \in \Theta\) such that \[ \hat{\theta} \in \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. To see what else it contains, expand \(\ell\) around \(\hat{\theta}\) to the second order: \[ \ell(\theta) \approx \ell(\hat{\theta}) - \frac{1}{2}\big(-\ell''(\hat{\theta})\big)(\theta - \hat{\theta})^2. \]
This is the quadratic approximation. 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:
- If \(-\ell''(\hat{\theta})\) is large, \(\ell\) drops off quickly, only \(\theta\) near \(\hat{\theta}\) explain the data well, and \(\hat{\theta}\) is precisely determined.
- If \(-\ell''(\hat{\theta})\) is small, \(\ell\) is flat, many values fit nearly as well, and \(\hat{\theta}\) is imprecise.
The curvature at \(\hat{\theta}\) is therefore exactly the right object to measure how informative the data are about \(\theta\). We will name it in a moment.
Example (normal with known variance). Let \(X_1, \ldots, X_n\) be iid \(N(\mu, \sigma^2)\) with \(\sigma^2\) known. The log-likelihood is \[ \begin{aligned} \ell(\mu) &= \sum_{i=1}^{n} \log\left[\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{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{aligned} \] By differentiating we get \(\ell'(\mu) = \frac{1}{\sigma^2}\sum_{i=1}^{n}(x_i - \mu)\), and setting \(\ell'(\mu) = 0\) gives the MLE \(\hat{\mu} = \bar{X}_n\), the sample mean.
Information and standard errors
The curvature at the MLE needs its own notation:
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 of \(I(\theta)\): the expected (or plain Fisher) information \[ \mathcal{I}(\theta) = E_\theta\big(I(\theta)\big) = \operatorname{Var}\big(U(\theta)\big), \] which averages the curvature over what we might have observed. The expected information is a property of the model and is computed once; the observed information is a property of this dataset and can be read off the log-likelihood we already have. For inference, we use the observed information at the MLE, because it reflects the curvature of the log-likelihood that we actually computed.
Standard error from curvature
We have now alluded many times to the role that information, as measured by the second derivative of the log-likelihood, plays when measuring precision of estimates — most recently in the innocent-looking final equation above. This follows from some relatively straightforward algebra, starting with \[ \begin{aligned} \int p_\theta(x)\,dx = 1 \quad\Longrightarrow\quad 0 &= \frac{d}{d\theta}\int p_\theta(x)\,dx = \int \frac{dp_\theta}{d\theta}\,dx \\ &\overset{(\dagger)}{=} \int p_\theta \frac{\partial \log p_\theta}{\partial\theta}\,dx = E_\theta\big(U(\theta)\big), \end{aligned} \tag{$\ast$} \] where \((\dagger)\) holds because \(\frac{dp}{d\theta} = p\,\frac{d\log p}{d\theta}\).
We can differentiate \(0 = E(U(\theta))\) once more to get \[ 0 = \frac{d}{d\theta}\int p_\theta U \,dx = \int \underbrace{\frac{dp_\theta}{d\theta}}_{= \,p_\theta U} U \,dx + \int p_\theta \underbrace{\frac{dU}{d\theta}}_{= \,-I} \,dx = E_\theta(U^2) - E\big(I(\theta)\big), \] so \(E_\theta(U^2) = E_\theta(I(\theta))\), and because \(E(U) = 0\) the left side is just \(\operatorname{Var}_\theta(U)\).
This argument works because we assume that we can switch the order of integration and differentiation — see Theorem 2.4.3 in CB for details — but this is true in most practical situations.
Next, consider the standard error of the MLE. We will begin by showing a nice heuristic for why it is reasonable to expect that the inverse of the observed information indeed is the variance of the MLE.
Compare the quadratic approximation \[ \ell(\theta) \approx \ell(\hat{\theta}) - \tfrac{1}{2} I(\hat{\theta})(\theta - \hat{\theta})^2 \] with the log-density of a \(N\big(\hat{\theta},\, I(\hat{\theta})^{-1}\big)\) distribution: \[ \log \varphi\big(\theta; \hat{\theta}, I(\hat{\theta})^{-1}\big) = \text{constant} - \tfrac{1}{2} I(\hat{\theta})(\theta - \hat{\theta})^2. \]
The two expressions agree up to an additive constant, meaning that the shape of the log-likelihood near the MLE looks exactly like a normal log-density with variance \(I(\hat{\theta})^{-1}\). Hence, it makes sense that inverting the curvature gives the variance and thus the standard deviation of the MLE: \[ SE(\hat{\theta}) \approx \frac{1}{\sqrt{I(\hat{\theta})}}. \qquad \blacksquare \]
A more precise argument for why this works. The quadratic-approximation argument above tells us that the likelihood looks normal-shaped near its peak with \(I(\hat{\theta})^{-1}\) as the variance of the approximating normal curve. But the standard error of \(\hat{\theta}\) is a statement about sampling variation. Why is the same number \(I(\hat{\theta})^{-1}\) the answer to both questions?
In the normal model with known variance the quadratic approximation is exact and there is nothing to check, but for a general regular model the link goes through a short asymptotic argument that we sketch below. It rests on two facts:
Fact 1. Variance of score \(=\) expected information, or, as stated and proved exactly above (for a single observation, \(n = 1\)): \[ \operatorname{Var}_\theta\big(U(\theta)\big) = E_\theta\big(I(\theta)\big) = \mathcal{I}(\theta). \]
Fact 2. The score satisfies a central limit theorem. Let \(\theta_0\) denote the true value of the parameter, the one that actually generated the data \(X_1, \ldots, X_n\). 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 random variables (one per observation). By Fact 1 applied with \(\theta = \theta_0\), each term has mean zero and variance \(\mathcal{I}(\theta_0)\) under \(p_{\theta_0}\), so by the CLT: \[ U(\theta_0)/\sqrt{n} \overset{d}{\to} N\big(0, \mathcal{I}(\theta_0)\big). \] It is the “mean zero” part that needs us to evaluate \(U\) at \(\theta = \theta_0\); under any other \(\theta\) (or \(p\)), the last equality of \((\ast)\) above would not hold.
Now, put the two facts together. The MLE solves \(U(\hat{\theta}) = 0\). Taylor-expanding \(U\) around the true value \(\theta_0\), \[ 0 = U(\hat{\theta}) \approx U(\theta_0) - I(\theta_0)(\hat{\theta} - \theta_0) \quad\Longrightarrow\quad \hat{\theta} - \theta_0 \approx \frac{U(\theta_0)}{I(\theta_0)}. \]
The numerator has variance approximately equal to \(n\mathcal{I}(\theta_0)\) (from Fact 2). The denominator is a sum of \(n\) iid terms with mean \(\mathcal{I}(\theta_0)\) (from Fact 1), so by the LLN it equals \(n\mathcal{I}(\theta_0)\) to leading order. Substituting, \[ \operatorname{Var}(\hat{\theta}) \approx \frac{n\mathcal{I}(\theta_0)}{\big(n\mathcal{I}(\theta_0)\big)^2} = \frac{1}{n\mathcal{I}(\theta_0)} \approx I(\theta_0)^{-1}, \] where we have been slightly sloppy with the notation, because we now refer to \(I(\theta_0)\) as the information contained in the whole sample, the expected value of which is \(n\) times the expected information \(\mathcal{I}(\theta_0)\) in a single observation.
Replacing \(\theta_0\) by \(\hat{\theta}\), which rests on the fundamental result that \(\hat{\theta}\) is a consistent estimator for \(\theta_0\) (not proven by us, yet), we arrive at the formula \[ SE(\hat{\theta}) = \frac{1}{\sqrt{I(\hat{\theta})}}. \]
Note that all “\(\approx\)”-signs can be replaced by proper asymptotic arguments in order to set up a rigorous proof of this result.
Example (normal with known variance, revisited). Continuing the example above, we have \[ \ell''(\mu) = -\frac{n}{\sigma^2}, \qquad I(\hat{\mu}) = -\ell''(\hat{\mu}) = \frac{n}{\sigma^2}. \] Therefore, \[ 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.
Several parameters
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, \] and 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}) = \underline{0}\), and the multivariate analogue of inverting the curvature gives \[ \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 the parameter estimates.
Example (normal distribution, both parameters unknown). 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 the 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 the partial derivatives to zero gives \[ \begin{aligned} \frac{\partial \ell}{\partial \mu} &= \frac{1}{\tau}\sum_{i=1}^{n}(x_i - \mu) = 0 &&\Longrightarrow\quad \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 &&\Longrightarrow\quad \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{X}_n)^2, \end{aligned} \] substituting \(\mu = \hat{\mu} = \bar{X}_n\) from the equation above.
Note the divisor in the MLE for \(\sigma^2\): it is \(n\), not \(n-1\). The MLE for \(\sigma^2\) is biased; maximum likelihood estimates are not always unbiased. (They are, however, asymptotically unbiased under some 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_{i=1}^{n}(x_i - \mu), \qquad \frac{\partial^2 \ell}{\partial\tau^2} = \frac{n}{2\tau^2} - \frac{1}{\tau^3}\sum_{i=1}^{n}(x_i - \mu)^2. \]
At the MLE, \(\sum_{i=1}^{n}(x_i - \mu) = 0\), which makes the off-diagonal elements vanish. Using \(\sum_i (x_i - \bar{X}_n)^2 = n\hat{\sigma}^2\), the observed information matrix evaluates to a diagonal matrix \[ 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, \[ SE(\hat{\mu}) \approx \frac{\hat{\sigma}}{\sqrt{n}}, \qquad SE(\hat{\sigma}^2) \approx \hat{\sigma}^2 \sqrt{2/n}. \]
The off-diagonal zero is not surprising, given the classical result that the sample mean and sample variance are statistically independent.
Example (linear regression with normal residuals). Let us apply the likelihood machinery to the ordinary simple linear regression, with residuals drawn from the normal distribution. A pleasant surprise (spoiler alert) is that the MLEs of the regression coefficients under normality are the same as the least squares estimates, which we can derive without assuming anything about the distribution of the residuals!
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{iid}{\sim} N(0, \sigma^2), \] with the \(x_i\) regarded as fixed. An equivalent way to write this model conditional on \(x_i\): \[ Y \sim N(\beta_0 + \beta_1 x_i,\, \sigma^2), \quad \text{independently.} \]
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 \[ SSR(\beta_0, \beta_1) = \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2. \]
In other words: the MLE of the regression coefficients is the ordinary least squares estimator. We obtain the familiar expressions by setting \(\partial\ell/\partial\beta_0 = 0\) and \(\partial\ell/\partial\beta_1 = 0\), leading to \[ \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, \] or, equivalently, \[ n\beta_0 + \beta_1 \sum_i x_i = \sum_i y_i, \qquad \beta_0 \sum_i x_i + \beta_1 \sum_i x_i^2 = \sum_i x_i y_i. \]
Solving this \(2 \times 2\) system gives the classical closed-form estimators. Writing \(\bar{x} = n^{-1}\sum_i x_i\), \(\bar{y} = n^{-1}\sum_i y_i\), \(S_{xx} = \sum_i (x_i - \bar{x})^2\), \(S_{xy} = \sum_i (x_i - \bar{x})(y_i - \bar{y})\), we get \[ \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} \big(y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i\big)^2 = \frac{SSR(\hat{\beta}_0, \hat{\beta}_1)}{n}. \]
The standard errors of the parameter estimates can be read off the inverse of the observed information matrix as before, but the algebra becomes heavier in this example.
Asymptotic normality of the MLE
We have already developed arguments that contain the substance of a limit theorem. To turn them into a formal result, we need to attach the correct \(\sqrt{n}\)-normalization and formulate the regularity conditions that justify each step. We also need the central limiting results from the previous lecture.
Theorem. 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 \(\theta\) 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 \overset{p}{\to} \theta_0\). Then \[ \sqrt{n}\big(\hat{\theta}_n - \theta_0\big) \overset{d}{\to} N\big(0, \mathcal{I}(\theta_0)^{-1}\big). \]
Sketch of proof. Apply Taylor’s Theorem to the score \(U\) around \(\theta_0\), in its mean-value form: 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). \]
We already know that \(U' = \ell'' = -I\), so we can rearrange and obtain \[ \hat{\theta}_n - \theta_0 = \frac{U(\theta_0)}{I(\theta_n^{*})}, \] and multiplying through by \(\sqrt{n}\), \[ \sqrt{n}\big(\hat{\theta}_n - \theta_0\big) = \frac{U(\theta_0)/\sqrt{n}}{I(\theta_n^{*})/n}. \]
We will treat the numerator and denominator separately and combine them afterwards with Slutsky’s Theorem.
Numerator. The sample-total score \(U(\theta_0) = \sum_i \frac{\partial}{\partial\theta}\log f_\theta(X_i)|_{\theta = \theta_0}\) is a sum of \(n\) iid random variables, each with mean \(0\) and variance \(\mathcal{I}(\theta_0)\). By the CLT that we proved in the previous lecture, \[ \frac{U(\theta_0)}{\sqrt{n}} \overset{d}{\to} N\big(0, \mathcal{I}(\theta_0)\big). \]
Denominator. The quantity \(I(\theta_n^{*})/n\) is a sample average of \(n\) iid terms, evaluated at the random point \(\theta_n^{*}\). Consistency of \(\hat{\theta}_n\) implies that \(\theta_n^{*} \overset{p}{\to} \theta_0\) (since \(\theta_n^{*}\) is squeezed between \(\hat{\theta}_n\) and \(\theta_0\)). Continuity (in \(\theta\)) of \(\frac{\partial^2}{\partial\theta^2}\log f_\theta(x)\) in a neighbourhood of \(\theta_0\), together with the law of large numbers, gives \[ \frac{I(\theta_n^{*})}{n} \overset{p}{\to} \mathcal{I}(\theta_0). \] Continuity is needed to control the behaviour of the random \(\theta_n^{*}\) relative to \(\theta_0\), but we skip the details.
Combine with Slutsky. Slutsky’s Theorem tells us that if the numerator converges in distribution and the denominator converges in probability to a positive constant, then the ratio converges in distribution: \[ \sqrt{n}\big(\hat{\theta}_n - \theta_0\big) = \frac{U(\theta_0)/\sqrt{n}}{I(\theta_n^{*})/n} \overset{d}{\to} \frac{N\big(0, \mathcal{I}(\theta_0)\big)}{\mathcal{I}(\theta_0)} = N\!\left(0, \frac{1}{\mathcal{I}(\theta_0)}\right). \] The last step uses the fact that dividing a \(N(0, \sigma^2)\) random variable by a positive constant \(a\) yields a \(N(0, \sigma^2/a^2)\) random variable. \(\qquad \blacksquare\)
A note on the consistency assumption. The proof above relied on the hypothesis that \(\hat{\theta}_n \overset{p}{\to} \theta_0\). Consistency is not automatic but follows from a fairly nice identity. By the LLN, the rescaled sample log-likelihood converges pointwise to a population analogue, \[ \frac{1}{n}\ell_n(\theta) \overset{p}{\to} M(\theta) = E_{\theta_0}\big[\log f_\theta(X)\big], \] and 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_{KL}(f_{\theta_0}, f_\theta) \geq 0, \] where this quantity is the Kullback–Leibler distance between \(f_{\theta_0}\) and \(f_\theta\), denoted \(D_{KL}(f_{\theta_0}, f_\theta)\), with equality if and only if \(\theta = \theta_0\) almost surely. The non-negativity is Jensen’s inequality applied to \(-\log(\cdot)\) (but we skip a few details here). The point is: under identifiability (distinct \(\theta\) giving distinct distributions), the limit \(M(\theta)\) is therefore strictly maximised at the true value \(\theta_0\).
So the expected log-likelihood is uniquely maximised at the truth, and the sample log-likelihood approximates it, and it logically follows (and it can be technically proven) that the sample maximizer \(\hat{\theta}_n\) must track the population maximizer \(\theta_0\).
The same identity also indicates what happens without a correctly specified model. If the true data generating distribution is not in the family \(\{f_\theta\}\), the KL inequality says the MLE still converges, but to the value \(\theta^{*}\) that minimizes \(D_{KL}(\text{truth}, f_\theta)\), so that the parameter, in this sense, makes the model the “closest” fit to a truth that does not belong to it. This can be the starting point of robustness analysis against model misspecification.
Curvature inversion, revisited. The theorem above can be restated as the approximation \[ \hat{\theta}_n \sim N\!\left(\theta_0, \frac{1}{n\mathcal{I}(\theta_0)}\right). \] Replacing \(n\mathcal{I}(\theta_0)\) by the observed information \(I(\theta_0)\) (by the LLN) and then \(\theta_0\) by \(\hat{\theta}_n\) (by consistency and continuous mapping), gives the working approximation \[ \hat{\theta}_n \sim N\!\left(\theta_0, I(\hat{\theta}_n)^{-1}\right). \]
This is the same statement that we used heuristically in the examples above to read standard errors off the curvature of \(\ell\), but 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\) we can use an entirely analogous argument (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)\)), which gives \[ \sqrt{n}\big(\hat{\theta}_n - \theta_0\big) \overset{d}{\to} N_p\big(\underline{0}, \mathcal{I}(\theta_0)^{-1}\big). \]
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}\big(g(\hat{\theta}_n) - g(\theta_0)\big) \overset{d}{\to} N\!\left(0, \frac{g'(\theta_0)^2}{\mathcal{I}(\theta_0)}\right), \] so that we immediately obtain the limit distribution of any smooth transformation of the MLE: odds ratios, log-odds, predicted probabilities, or anything else that we might want to report.
The following result is not very surprising in that respect:
Theorem 7.2.10 in CB. If \(\hat{\theta}\) is the MLE of \(\theta\), then \(\tau(\hat{\theta})\) is the MLE of \(\tau(\theta)\).