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: the observations come from a model that is known up to a set of parameters, such as \(X_1, \ldots, X_n \overset{\text{iid}}{\sim} N(\mu, \sigma^2)\), or \(Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\) with \(\varepsilon_i \sim N(0, \sigma^2)\). The task is to estimate the parameters using data. Estimating \(f(\cdot)\) directly, without such assumptions, is called nonparametric statistics.
We also concentrate on frequentist statistics, where parameters are fixed but unknown numbers, and the data are represented by random variables. In Bayesian statistics it is, in a sense, the other way around.
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, and that \(E(\bar{X}) = \mu\). The sample mean is therefore a very natural estimator of \(\mu\).
2. The method of moments. 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 solve the system \[m_k = \mu_k(\theta_1, \ldots, \theta_p), \quad k = 1, \ldots, p.\] For a normal sample, this gives \(\tilde{\theta} = \bar{X}\) and \(\tilde{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2\) (Example 7.2.1, see the notes).
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\).
The likelihood is not a probability distribution over \(\Theta\) and is not normalised. But the ratio \(L(\theta_2)/L(\theta_1)\) tells us how many times more probable the data are under \(\theta_2\) than under \(\theta_1\). For a continuous model we use the density, \(L(\theta;\, x) = p_\theta(x)\), which is the same ranking up to a constant that does not involve \(\theta\) (see the notes).
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).\]
Once we have \(\ell(\theta)\), the data are essentially summarised. The point estimate (where is \(\ell\) maximised?) and the standard error (how sharply is \(\ell\) peaked there?) 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.
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.
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\).
The coin example with \(n\) tosses and a proportion \(\hat\theta\) of heads. Left: the likelihood relative to its maximum. Right: the log-likelihood and its quadratic approximation (dashed). The readout shows the information \(I(\hat\theta) = -\ell''(\hat\theta)\), defined on the next slide.
Definition (Score and observed information). The score function is \(U(\theta) = \ell'(\theta)\), and the observed Fisher information is \[I(\theta) = -\ell''(\theta).\]
This suggests measuring the precision of \(\hat\theta\) by the following formula, where more information means smaller standard error: \[\operatorname{SE}(\hat\theta) \approx \frac{1}{\sqrt{I(\hat\theta)}}.\]
For an iid sample, the log-likelihood is a sum, \(\ell(\theta) = \sum_{i=1}^n \log p_\theta(X_i)\), so its curvature is a sum of \(n\) contributions, one from each observation.
The expected Fisher information \(\mathcal{I}(\theta)\) is the average contribution of one observation: \[\mathcal{I}(\theta) = E_\theta\!\left(-\frac{\partial^2}{\partial\theta^2}\log p_\theta(X)\right), \qquad E_\theta\bigl(I(\theta)\bigr) = n\,\mathcal{I}(\theta).\]
Hence, if we believe the formula from the previous slide (will be formalized later): \[\operatorname{SE}(\hat\theta) \approx \frac{1}{\sqrt{I(\hat\theta)}} \approx \frac{1}{\sqrt{n\,\mathcal{I}(\theta)}},\] which is the \(\sqrt{n}\) rate that we recognize from the CLT.
| Model | \(\mathcal{I}(\theta)\), one observation | \(\operatorname{SE}\) of the MLE |
|---|---|---|
| Bernoulli\((\theta)\) | \(\dfrac{1}{\theta(1-\theta)}\) | \(\sqrt{\hat\theta(1-\hat\theta)/n}\) |
| \(N(\mu, \sigma^2)\), \(\sigma\) known | \(\dfrac{1}{\sigma^2}\) | \(\sigma/\sqrt{n}\) |
| Exponential, rate \(\lambda\) | \(\dfrac{1}{\lambda^2}\) | \(\hat\lambda/\sqrt{n}\) |
Theorem (Asymptotic normality of the MLE). Under regularity conditions (see the notes), \[\sqrt{n}\,(\hat\theta_n - \theta_0) \xrightarrow{d} N\!\bigl(0,\, \mathcal{I}(\theta_0)^{-1}\bigr),\] where \(\theta_0\) is the true parameter value.
In the language of limit theory: the asymptotic variance of the MLE is \(1/\mathcal{I}(\theta_0)\), and \(\hat\theta_n\) is approximately \(N\big(\theta_0,\ 1/(n\mathcal{I}(\theta_0))\big)\).
The idea of the proof (more details in the notes):
The asymptotic variance depends on the unknown \(\theta_0\). Replacing \(n\mathcal{I}(\theta_0)\) by the observed information \(I(\hat\theta_n)\) gives \[\hat\theta_n \;\stackrel{\textrm{appr}}{\sim}\; N\!\bigl(\theta_0,\; I(\hat\theta_n)^{-1}\bigr), \qquad \operatorname{SE}(\hat\theta_n) \approx \frac{1}{\sqrt{I(\hat\theta_n)}},\] which is the “idea” from earlier about what is a sensible measure of uncertainty of the MLE.
We can then construct an approximate (i.e., correct in the limit) confidence interval \[\hat\theta_n \pm z_{\alpha/2}\,\frac{1}{\sqrt{I(\hat\theta_n)}}\] with coverage tending to \(1 - \alpha\) as \(n\rightarrow \infty\). It is the workhorse of large-sample inference.
Let \(X_1, \ldots, X_n\) be iid exponential with rate \(\lambda\), so \(p_\lambda(x) = \lambda e^{-\lambda x}\) for \(x > 0\). Then \[\ell(\lambda) = n\log\lambda - \lambda\sum_{i=1}^n X_i, \qquad \ell'(\lambda) = \frac{n}{\lambda} - \sum_{i=1}^n X_i, \qquad I(\lambda) = -\ell''(\lambda) = \frac{n}{\lambda^2}.\]
The MLE is \(\hat\lambda = 1/\bar{X}_n\), and \[\operatorname{SE}(\hat\lambda) \approx \frac{1}{\sqrt{I(\hat\lambda)}} = \frac{\hat\lambda}{\sqrt{n}}.\] Here \(\mathcal{I}(\lambda) = 1/\lambda^2\), so the asymptotic variance is \(\lambda^2\): the same as we found with the delta method in Topic 5.
With \(n = 50\) times between customer arrivals and \(\bar{x} = 2.5\) minutes: \(\hat\lambda = 0.4\) per minute, \(\operatorname{SE} = 0.4/\sqrt{50} = 0.057\), and a 95% interval for \(\lambda\) is \(0.4 \pm 1.96 \cdot 0.057 = [0.29,\ 0.51]\).
Illustration of the main result. Each repetition draws \(n\) observations from the exponential distribution with rate \(\lambda_0 = 0.4\) and computes \(\hat\lambda = 1/\bar{X}_n\). Left: 4000 values of \(\sqrt{n}(\hat\lambda - \lambda_0)\) against the \(N(0, \lambda_0^2)\) density. Right: 50 of the intervals \(\hat\lambda \pm 1.96\,\hat\lambda/\sqrt{n}\).
For a parameter vector \(\theta = (\theta_1, \ldots, \theta_p)\), the score is the gradient \(U(\theta) = \nabla\ell(\theta)\), and the observed information matrix is the negative Hessian \(I(\theta) = -\nabla^2\ell(\theta)\). The MLE solves \(U(\hat\theta) = \mathbf{0}\), and \[\widehat{\operatorname{Cov}}(\hat\theta) \approx I(\hat\theta)^{-1}.\] The diagonal gives the variances of the individual estimates, and the off-diagonal elements their covariances.
Example: normal, both parameters unknown. The MLEs and the inverse information matrix are \[\hat\mu = \bar{X}_n, \qquad \hat\sigma^2 = \frac{1}{n}\sum_{i=1}^n(X_i - \bar{X}_n)^2, \qquad I(\hat\theta)^{-1} = \begin{pmatrix} \hat\sigma^2/n & 0 \\ 0 & 2\hat\sigma^4/n \end{pmatrix}.\]
The derivation is in the notes.
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\).
A pleasant surprise: the least squares estimates can be derived without assuming anything about the distribution of the residuals!
Theorem 7.2.10 (Invariance). If \(\hat\theta\) is the MLE of \(\theta\), then \(\tau(\hat\theta)\) is the MLE of \(\tau(\theta)\).
We do not have to re-derive the MLE for a transformed parameter, and the delta method from Topic 5 gives its asymptotic distribution: if \(g'(\theta_0) \neq 0\), \[\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).\]
This gives the point estimate and the standard error of any smooth function of the parameters (odds ratios, log-odds, predicted probabilities), directly from the MLE and the observed information. With several parameters, the same holds with the information matrix and the multivariate delta method.
| Question | Answer from the likelihood |
|---|---|
| Which \(\theta\) fits the data best? | The MLE, where \(\ell\) peaks: \(U(\hat\theta) = 0\) |
| How precise is it? | The information, the curvature at the peak: \(\operatorname{SE}(\hat\theta) \approx 1/\sqrt{I(\hat\theta)}\) |
| How does precision grow with \(n\)? | Information adds up, \(I \approx n\,\mathcal{I}(\theta)\), so the SE shrinks like \(1/\sqrt{n}\) |
| Why does that work? | \(\sqrt{n}(\hat\theta_n - \theta_0) \xrightarrow{d} N(0, \mathcal{I}(\theta_0)^{-1})\): the CLT and the LLN applied to the log-likelihood (proof in the notes) |
| What do we report? | \(\hat\theta \pm z_{\alpha/2}/\sqrt{I(\hat\theta)}\), and the delta method for functions of \(\theta\) |