Lesson overview | Lesson overview | Next part

Estimation Theory: Part 1: Intuition to 5. Maximum Likelihood Estimation

1. Intuition

1.1 What Is Statistical Estimation?

Suppose you flip a coin 100 times and observe 63 heads. The coin's true bias $p$ - the probability of heads on any single flip - is unknown. What is your best guess for $p$? The obvious answer is $\hat{p} = 0.63$. But why is this the right answer? What makes it "best"? Could a different estimate be better in some sense? And how confident should you be - could the true $p$ be 0.5 (a fair coin) despite observing 63 heads?

These are precisely the questions estimation theory addresses. The statistical estimation problem has three ingredients:

1. An unknown parameter $\theta \in \Theta$ governing some probability distribution $p(\mathbf{x}; \theta)$
2. Observed data $\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}$ drawn independently from $p(\cdot; \theta)$
3. An estimator $\hat{\theta} = T(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)})$: a function of the data that produces a guess for $\theta$

The key insight, often overlooked by beginners, is that the estimator is a random variable. Before we collect data, $\hat{\theta}$ is uncertain - its value depends on which random sample we happen to observe. If we repeated the experiment with a fresh sample, we would get a different $\hat{\theta}$. The sampling distribution of $\hat{\theta}$ - the distribution of $\hat{\theta}$ across all possible samples - is the object of study. A good estimator has a sampling distribution tightly concentrated around the true $\theta$.

For AI: This perspective matters directly. When you train a neural network on a dataset $\mathcal{D}$ and obtain weights $\hat{\boldsymbol{\theta}}$, those weights are a realisation of an estimator - a function of the random training data. Training on a different shuffle or subsample gives different $\hat{\boldsymbol{\theta}}$. The variation you see across training runs, seeds, and dataset splits is the sampling distribution of the MLE manifesting in practice.

1.2 The Estimation Landscape

Estimation theory spans three major frameworks, each answering a different question:

Point estimation produces a single value $\hat{\theta}$ for the unknown parameter. The challenge is specifying what "best" means. Three criteria dominate:

• Unbiasedness: $\mathbb{E}[\hat{\theta}] = \theta$ - on average, the estimator is correct
• Minimum variance: among unbiased estimators, prefer the one with smallest $\operatorname{Var}(\hat{\theta})$
• Consistency: $\hat{\theta}_n \xrightarrow{P} \theta$ as $n \to \infty$ - with enough data, the estimator converges to the truth

These criteria can conflict. The Cramer-Rao lower bound establishes a hard floor on variance for unbiased estimators, and the estimators that achieve this floor are called efficient. Maximum likelihood estimation (MLE) is the dominant method: it is consistent, asymptotically efficient, and invariant under smooth reparametrisations.

Interval estimation goes beyond a point to a range $[\hat{\theta}_L, \hat{\theta}_U]$ with a coverage guarantee: in repeated experiments, $1 - \alpha$ fraction of such intervals will contain the true $\theta$. Frequentist confidence intervals capture estimation uncertainty without requiring a prior distribution on $\theta$.

Asymptotic theory studies estimator behaviour as $n \to \infty$. The flagship result is the asymptotic normality of MLE: under regularity conditions,

$$\sqrt{n}\left(\hat{\theta}_{\text{MLE}} - \theta^*\right) \xrightarrow{d} \mathcal{N}\!\left(0,\, \mathcal{I}(\theta^*)^{-1}\right)$$

where $\mathcal{I}(\theta^*)$ is the Fisher information. This single result simultaneously justifies MLE as an estimation method, provides the basis for asymptotic confidence intervals, and connects estimation theory to information geometry.

For AI: The three frameworks map directly to ML practice. Point estimation = model training (find $\hat{\boldsymbol{\theta}}$). Interval estimation = evaluation with error bars (report accuracy $\pm$ 95% CI). Asymptotic theory = understanding why larger datasets produce more reliable models (the $1/\sqrt{n}$ convergence rate).

1.3 Historical Timeline

ESTIMATION THEORY - HISTORICAL TIMELINE
========================================================================

  1795  Gauss (age 18) uses least squares to predict asteroid orbits
        - the first systematic estimation method

  1805  Legendre publishes least squares formally (Gauss priority dispute)

  1809  Gauss derives least squares from the assumption of Gaussian errors
        - first connection between MLE and squared-error minimisation

  1894  Pearson introduces method of moments - systematic moment matching

  1912  Fisher introduces "maximum likelihood" as a term; begins systematic
        study of estimation properties

  1922  Fisher: "On the Mathematical Foundations of Theoretical Statistics"
        - defines sufficiency, efficiency, consistency; derives properties of MLE

  1925  Fisher introduces Fisher information and the information inequality
        (precursor to Cramer-Rao)

  1945  Cramer proves the Cramer-Rao lower bound (independently of Rao)

  1945  Rao independently proves the same bound + Rao-Blackwell theorem

  1947  Lehmann & Scheffe: completeness and UMVUE theory

  1950s Wald develops statistical decision theory - unified framework for
        estimation as minimising expected loss

  1979  Efron introduces the bootstrap - non-parametric confidence intervals
        without distributional assumptions

  1998  Amari: natural gradient via Fisher information geometry
        - direct connection to modern ML optimisation

  2015  Martens & Grosse: K-FAC - tractable approximation of FIM for DNNs

  2017  Kirkpatrick et al.: Elastic Weight Consolidation (EWC) uses FIM
        to prevent catastrophic forgetting in continual learning

  2022  Hoffmann et al. (Chinchilla): scaling laws estimated via MLE on
        (N, D, L) data - estimation theory at trillion-parameter scale

========================================================================

1.4 Why Estimation Theory Matters for AI

Every major component of a modern ML pipeline rests on estimation theory:

Training objectives. Minimising cross-entropy loss is equivalent to performing MLE. The connection is exact: $\arg\min_{\boldsymbol{\theta}} \sum_i -\log p(y^{(i)} \mid \mathbf{x}^{(i)}; \boldsymbol{\theta}) = \arg\max_{\boldsymbol{\theta}} \prod_i p(y^{(i)} \mid \mathbf{x}^{(i)}; \boldsymbol{\theta})$. Every training run of GPT-4, LLaMA-3, Gemini, or any neural network trained with NLL loss is performing MLE on the conditional distribution.

Evaluation uncertainty. When a benchmark reports "Model A achieves 73.2% accuracy", this is a point estimate. How uncertain is it? With $n = 1000$ test examples, the 95% CI is approximately $\pm 2.8\%$ - meaning Model A and Model B with 73.5% accuracy are statistically indistinguishable. Estimation theory makes this precise.

Second-order optimisation. Adam, K-FAC, and natural gradient all incorporate curvature information. The natural gradient multiplies the gradient by $\mathcal{I}(\boldsymbol{\theta})^{-1}$, where $\mathcal{I}$ is the Fisher information matrix. This reparametrises the optimisation in the geometry of the statistical manifold, making gradient steps invariant to reparametrisation. K-FAC (Martens & Grosse, 2015) approximates $\mathcal{I}^{-1}$ tractably for deep networks.

Continual learning. Elastic weight consolidation (Kirkpatrick et al., 2017) protects important weights during sequential task learning. "Importance" is measured by the diagonal of the FIM: parameters with high Fisher information are those to which the likelihood is most sensitive, and perturbing them most damages performance.

Calibration. A well-calibrated model's confidence scores match empirical frequencies: when it says "80% confident", it's right 80% of the time. Temperature scaling - dividing logits by a scalar $\tau$ - is an MLE problem: find the $\tau$ that maximises likelihood on a held-out calibration set.

---

2. The Formal Estimation Problem

2.1 Statistical Model and Parametric Family

Definition 2.1 (Parametric Statistical Model). A parametric statistical model is a collection of probability distributions indexed by a parameter:

$$\mathcal{M} = \{p(\cdot\,;\boldsymbol{\theta}) : \boldsymbol{\theta} \in \Theta\}$$

where $\Theta \subseteq \mathbb{R}^k$ is the parameter space. The model is:

• Correctly specified if the true data-generating distribution $p^*$ satisfies $p^* = p(\cdot\,;\boldsymbol{\theta}^*)$ for some $\boldsymbol{\theta}^* \in \Theta$
• Misspecified if $p^* \notin \mathcal{M}$ - the model family does not contain the truth (most practical ML models are misspecified)
• Identifiable if different parameters give different distributions: $\boldsymbol{\theta} \neq \boldsymbol{\theta}' \Rightarrow p(\cdot\,;\boldsymbol{\theta}) \neq p(\cdot\,;\boldsymbol{\theta}')$

Identifiability is necessary for consistent estimation - you cannot consistently estimate a parameter that leaves the data distribution unchanged.

Standard examples of parametric families:

Non-examples (non-identifiable):

• $\mathcal{N}(\mu_1 + \mu_2, 1)$ with parameter $(\mu_1, \mu_2)$: we can only estimate $\mu_1 + \mu_2$, not the individual components. The model is not identifiable.
• A neural network with two hidden layers of width 1 using $\tanh$: swapping the two neurons gives the same function, so the parameterisation is not identifiable (though the function class may be).

The iid assumption. Throughout this section, observations $\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}$ are assumed independent and identically distributed (iid) from $p(\cdot\,;\boldsymbol{\theta}^*)$. The iid assumption enables the joint log-likelihood to factorise as a sum:

$$\log p(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)};\boldsymbol{\theta}) = \sum_{i=1}^n \log p(\mathbf{x}^{(i)};\boldsymbol{\theta})$$

This factorisation is what makes MLE computationally tractable and theoretically tractable. In practice, this assumption is violated whenever data points are correlated (time series, text sequences, spatially clustered data), and estimation theory has extensions to handle dependent data.

2.2 Point Estimators

Definition 2.2 (Estimator). An estimator of $\boldsymbol{\theta}$ is any measurable function $\hat{\boldsymbol{\theta}}_n = T(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)})$ of the sample. The key points:

1. $\hat{\boldsymbol{\theta}}_n$ is a random variable - before observing data, it is uncertain
2. A specific value $T(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}) = t$ after observing data is called an estimate (not estimator)
3. The sampling distribution of $\hat{\boldsymbol{\theta}}_n$ is the distribution of $T$ across all possible samples of size $n$

This distinction - estimator (random variable) vs. estimate (specific realisation) - is pedantic but consequential. Confidence interval coverage is a property of the estimator (random), not the estimate (fixed). When we say "95% CI", we mean that the random interval contains $\boldsymbol{\theta}^*$ with probability 95%, not that this specific computed interval has a 95% chance of containing $\boldsymbol{\theta}^*$ (after observing data, $\boldsymbol{\theta}^*$ is either in the interval or not - there is no remaining randomness).

Examples of estimators for the Gaussian mean $\mu$:

Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} \mathcal{N}(\mu, \sigma^2)$.

• Sample mean: $\hat{\mu}_1 = \frac{1}{n}\sum_{i=1}^n x^{(i)}$ - unbiased, consistent, efficient
• Constant zero: $\hat{\mu}_2 = 0$ - biased unless $\mu = 0$, inconsistent
• First observation: $\hat{\mu}_3 = x^{(1)}$ - unbiased! but inconsistent (variance $\sigma^2$ regardless of $n$)
• Trimmed mean: $\hat{\mu}_4 =$ mean of middle 90% - slightly biased for Gaussian, but robust to outliers
• Constant $c$: $\hat{\mu}_5 = c$ - biased by $|c - \mu|$, inconsistent; can have lower MSE than the sample mean when $|\mu - c|$ is small (illustrates bias-variance trade-off)

The existence of unbiased but inconsistent estimators ($\hat{\mu}_3$) and consistent but biased estimators (trimmed mean) shows these properties are logically independent.

2.3 The Loss Framework: MSE, Bias, and Variance

To compare estimators, we need a criterion. The most widely used is Mean Squared Error (MSE):

Definition 2.3 (MSE, Bias, Variance). For an estimator $\hat{\theta}$ of a scalar parameter $\theta$:

$$\operatorname{MSE}(\hat{\theta}) = \mathbb{E}\left[(\hat{\theta} - \theta)^2\right]$$

The bias-variance decomposition is the central identity of estimation theory:

$$\operatorname{MSE}(\hat{\theta}) = \operatorname{Bias}^2(\hat{\theta}) + \operatorname{Var}(\hat{\theta})$$

Proof:

$$\operatorname{MSE}(\hat{\theta}) = \mathbb{E}\left[(\hat{\theta} - \theta)^2\right]$$

Add and subtract $\mathbb{E}[\hat{\theta}]$:

$$= \mathbb{E}\left[\left((\hat{\theta} - \mathbb{E}[\hat{\theta}]) + (\mathbb{E}[\hat{\theta}] - \theta)\right)^2\right]$$ $$= \mathbb{E}\left[(\hat{\theta} - \mathbb{E}[\hat{\theta}])^2\right] + 2(\mathbb{E}[\hat{\theta}] - \theta)\underbrace{\mathbb{E}[\hat{\theta} - \mathbb{E}[\hat{\theta}]]}_{=0} + (\mathbb{E}[\hat{\theta}] - \theta)^2$$ $$= \operatorname{Var}(\hat{\theta}) + \operatorname{Bias}(\hat{\theta})^2 \qquad \square$$

where $\operatorname{Bias}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta$.

Geometric interpretation:

BIAS-VARIANCE DECOMPOSITION - GEOMETRIC VIEW
========================================================================

  True parameter \\theta* = 0 (target)

  Case 1: Low bias, low variance (good estimator)
    ---    Estimates cluster tightly around \\theta*
    ----   MSE \\approx small
    ---

  Case 2: Low bias, high variance (unbiased but noisy)
        -             -
    -       -   x   -       -
        -             -
    Centred on \\theta* but spread out.  MSE = Var is large.

  Case 3: High bias, low variance (biased but stable)
              ---
             ----    x  (\\theta*)
              ---
    Centred off target.  MSE = Bias^2 + small Var.

  Case 4: High bias, high variance (worst case)
    -           -
            -       x  (\\theta*)
        -       -
    Both centred off target AND spread out.

  Key trade-off: Reducing variance by shrinkage introduces bias.
  Increasing bias via regularisation often reduces variance more than
  it increases bias-squared, giving lower MSE overall.

========================================================================

The bias-variance trade-off in ML is the exact same phenomenon: regularisation (L2, dropout, early stopping) introduces bias (the model is pulled away from the MLE toward a constrained region) but reduces variance (the model is less sensitive to the particular training sample), often reducing test MSE overall.

Minimax estimators. MSE is not the only loss criterion. The minimax estimator minimises the worst-case risk: $\min_T \max_\theta \mathbb{E}[(\hat{\theta} - \theta)^2]$. The James-Stein estimator (1961) famously showed that in $\mathbb{R}^k$ for $k \geq 3$, the sample mean is inadmissible - there exists an estimator with strictly lower MSE at every $\theta$. This is the Stein paradox: shrinkage toward the origin (a biased estimator!) uniformly dominates the unbiased sample mean in 3+ dimensions.

2.4 Examples of Biased and Unbiased Estimators

Example 1: Sample mean is unbiased.

Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} p(\cdot;\theta)$ with $\mathbb{E}[x^{(i)}] = \mu$. Then:

$$\mathbb{E}\!\left[\bar{x}\right] = \mathbb{E}\!\left[\frac{1}{n}\sum_{i=1}^n x^{(i)}\right] = \frac{1}{n}\sum_{i=1}^n \mathbb{E}[x^{(i)}] = \frac{n\mu}{n} = \mu$$

The sample mean is an unbiased estimator of $\mu$ for any distribution with finite mean.

Example 2: MLE of Gaussian variance is biased.

Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} \mathcal{N}(\mu, \sigma^2)$. The MLE of $\sigma^2$ is:

$$\hat{\sigma}^2_{\text{MLE}} = \frac{1}{n}\sum_{i=1}^n (x^{(i)} - \bar{x})^2$$

Computing the bias: Using $x^{(i)} - \bar{x} = (x^{(i)} - \mu) - (\bar{x} - \mu)$:

$$\sum_{i=1}^n (x^{(i)} - \bar{x})^2 = \sum_{i=1}^n (x^{(i)} - \mu)^2 - n(\bar{x} - \mu)^2$$

Taking expectations:

$$\mathbb{E}\!\left[\sum_{i=1}^n (x^{(i)} - \bar{x})^2\right] = n\sigma^2 - n \cdot \frac{\sigma^2}{n} = (n-1)\sigma^2$$

Therefore:

$$\mathbb{E}[\hat{\sigma}^2_{\text{MLE}}] = \frac{(n-1)\sigma^2}{n} = \sigma^2 - \frac{\sigma^2}{n}$$

The bias is $-\sigma^2/n$ - the MLE underestimates variance. The correction is Bessel's correction: $s^2 = \frac{1}{n-1}\sum(x^{(i)} - \bar{x})^2$ is unbiased.

Recall from Section01: We saw Bessel's correction $s^2 = \frac{1}{n-1}\sum(x^{(i)}-\bar{x})^2$ introduced as the "standard" sample variance formula. Now we understand why: MLE gives $1/n$, which is biased. The $n-1$ denominator corrects for the one degree of freedom lost in estimating $\mu$ by $\bar{x}$.

Example 3: MLE is biased but consistent.

The Gaussian variance MLE $\hat{\sigma}^2_{\text{MLE}}$ has bias $-\sigma^2/n \to 0$ as $n \to \infty$. It is asymptotically unbiased and consistent. This illustrates that bias can be acceptable if it vanishes with sample size.

Example 4: When bias reduces MSE.

Consider estimating $\mu$ with the shrinkage estimator $\hat{\mu}_c = c \bar{x}$ for some $c \in (0,1)$. Then:

$$\operatorname{MSE}(\hat{\mu}_c) = (c-1)^2 \mu^2 + c^2 \frac{\sigma^2}{n}$$

The optimal $c^* = \frac{\mu^2}{\mu^2 + \sigma^2/n}$, which gives lower MSE than $c = 1$ (unbiased sample mean) whenever $\sigma^2/n > 0$. The key insight: when $|\mu|$ is small relative to $\sigma/\sqrt{n}$, shrinking toward zero introduces small bias but large variance reduction.

---

3. Properties of Estimators

3.1 Bias and Unbiasedness

Definition 3.1 (Bias). The bias of an estimator $\hat{\theta}$ for parameter $\theta$ is:

$$\operatorname{Bias}(\hat{\theta};\theta) = \mathbb{E}[\hat{\theta}] - \theta$$

The estimator is unbiased if $\operatorname{Bias}(\hat{\theta};\theta) = 0$ for all $\theta \in \Theta$, and asymptotically unbiased if $\operatorname{Bias}(\hat{\theta}_n;\theta) \to 0$ as $n \to \infty$.

Why unbiasedness is desirable - but not sacred. An unbiased estimator is correct on average over infinite repetitions of the experiment. This is a natural property to want. However, unbiasedness alone is insufficient:

• The estimator $\hat{\theta} = x^{(1)}$ (use only the first observation) is unbiased but throws away $n-1$ observations - clearly wasteful
• No unbiased estimator exists for some parameters: there is no unbiased estimator of $e^\mu$ based on a Gaussian sample (by the Lehmann-Scheffe theorem, for this to exist $e^\mu$ would need to be a polynomial in $\bar{x}$, which it is not)
• The Stein phenomenon (Section2.3) shows that in 3+ dimensions, unbiased estimators can be uniformly dominated by biased ones

Bias correction. When an estimator $\hat{\theta}$ has known bias $b(\theta) = \mathbb{E}[\hat{\theta}] - \theta$, we can correct it: $\hat{\theta}_{\text{corr}} = \hat{\theta} - b(\hat{\theta})$. This is the idea behind Bessel's correction ($n-1$ vs $n$ in the variance formula) and jackknife bias reduction.

For AI: Batch normalisation computes $\hat{\mu}_B = \frac{1}{B}\sum_{i \in \text{batch}} x_i$ as the batch mean during training. This is an unbiased estimator of the feature mean. During inference, PyTorch uses the running mean accumulated over training, which is a biased estimator of the current data mean if the distribution has shifted - illustrating that unbiasedness depends on whether the estimation scenario matches the training scenario.

3.2 Consistency

Definition 3.2 (Consistency). An estimator $\hat{\theta}_n$ is:

• Weakly consistent if $\hat{\theta}_n \xrightarrow{P} \theta$ for all $\theta \in \Theta$ (convergence in probability)
• Strongly consistent if $\hat{\theta}_n \xrightarrow{a.s.} \theta$ for all $\theta \in \Theta$ (almost-sure convergence)

Weak consistency: for every $\varepsilon > 0$, $P(|\hat{\theta}_n - \theta| > \varepsilon) \to 0$ as $n \to \infty$.

Sufficient conditions for consistency:

• If $\operatorname{Bias}(\hat{\theta}_n;\theta) \to 0$ and $\operatorname{Var}(\hat{\theta}_n) \to 0$ as $n \to \infty$, then $\hat{\theta}_n$ is weakly consistent. (By Chebyshev + the bias decomposition: $P(|\hat{\theta}_n - \theta| > \varepsilon) \leq \operatorname{MSE}(\hat{\theta}_n)/\varepsilon^2 \to 0$.)

Example: Sample mean consistency.

$\bar{x}_n = \frac{1}{n}\sum_{i=1}^n x^{(i)}$ satisfies $\operatorname{Bias} = 0$ and $\operatorname{Var}(\bar{x}_n) = \sigma^2/n \to 0$. Therefore $\bar{x}_n \xrightarrow{P} \mu$. This is exactly the Weak Law of Large Numbers - consistency of the sample mean is the LLN.

Example: MLE of Gaussian variance is consistent.

$\hat{\sigma}^2_{\text{MLE}} = \frac{1}{n}\sum(x^{(i)} - \bar{x})^2$ has bias $-\sigma^2/n \to 0$ and variance $O(1/n)$, so it is consistent even though it is biased for finite $n$.

Inconsistent estimator example: $\hat{\theta} = x^{(1)}$ (take only the first observation) has $\operatorname{Var} = \sigma^2$ regardless of $n$. It never concentrates - it is inconsistent.

Why consistency matters: Consistency is the minimal requirement for an estimator to be scientifically useful. A method that doesn't converge to the right answer with infinite data is fundamentally broken. Consistency does not guarantee that the estimator is good for small $n$ - convergence can be arbitrarily slow - but it at least ensures the method is correct in the large-data limit.

For AI: The "double descent" phenomenon in over-parameterised neural networks challenges classical bias-variance thinking. A massively over-parameterised model ($d \gg n$) can still be consistent for the Bayes-optimal predictor if trained with appropriate implicit regularisation (e.g., gradient descent from zero initialisation). This is an active area connecting statistical estimation theory to modern deep learning theory.

3.3 Efficiency

Among all unbiased estimators of $\theta$, which has the smallest variance? This question is answered by the Cramer-Rao bound (Section 4), but we can define efficiency as a property here.

Definition 3.3 (Efficiency and Relative Efficiency).

• An unbiased estimator $\hat{\theta}$ is efficient if its variance equals the Cramer-Rao lower bound: $\operatorname{Var}(\hat{\theta}) = 1/\mathcal{I}(\theta)$.
• The relative efficiency of two unbiased estimators $\hat{\theta}_1$ and $\hat{\theta}_2$ is $e(\hat{\theta}_1, \hat{\theta}_2) = \operatorname{Var}(\hat{\theta}_2)/\operatorname{Var}(\hat{\theta}_1)$. Values greater than 1 mean $\hat{\theta}_1$ is more efficient.

Examples:

1. Gaussian mean: The sample mean $\bar{x}$ is efficient - it achieves the CRB $\sigma^2/n$. Relative efficiency of the sample median vs. mean is $2/\pi \approx 0.637$ for Gaussian data - the median throws away ~36% of efficiency.

2. Gaussian variance: The unbiased sample variance $s^2$ is not efficient (the MLE $\hat{\sigma}^2$ achieves the CRB but is biased; correcting for bias loses the CRB).

Asymptotic efficiency. For large $n$, MLE achieves the CRB asymptotically - it is asymptotically efficient. This is one of its key advantages.

For AI: The sample mean is the most efficient estimator of the mean for Gaussian data. Batch gradient descent using all $n$ samples computes the exact gradient (equivalent to efficient estimation), while SGD uses a single sample or mini-batch (less efficient, higher variance, but much faster per update). The trade-off between statistical efficiency and computational efficiency is fundamental to modern ML training.

3.4 Sufficiency and the Rao-Blackwell Theorem

A sufficient statistic captures all the information in the data that is relevant to estimating $\theta$.

Definition 3.4 (Sufficient Statistic). A statistic $T = T(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)})$ is sufficient for $\theta$ if the conditional distribution $p(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)} \mid T = t; \theta)$ does not depend on $\theta$ for any $t$.

Intuitively: once you know $T$, the raw data $\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}$ provide no additional information about $\theta$.

Fisher-Neyman Factorisation Theorem. $T$ is sufficient for $\theta$ if and only if the likelihood factors as:

$$p(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)};\theta) = g(T(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}), \theta) \cdot h(\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)})$$

where $g$ depends on data only through $T$, and $h$ does not depend on $\theta$.

Examples of sufficient statistics:

Rao-Blackwell Theorem. If $\hat{\theta}$ is any unbiased estimator and $T$ is a sufficient statistic, then $\tilde{\theta} = \mathbb{E}[\hat{\theta} \mid T]$ satisfies:

$$\operatorname{Var}(\tilde{\theta}) \leq \operatorname{Var}(\hat{\theta})$$

with equality iff $\hat{\theta}$ is already a function of $T$. The "Rao-Blackwellised" estimator is at least as good.

Proof sketch: By the law of total variance, $\operatorname{Var}(\hat{\theta}) = \operatorname{Var}(\mathbb{E}[\hat{\theta}|T]) + \mathbb{E}[\operatorname{Var}(\hat{\theta}|T)] \geq \operatorname{Var}(\tilde{\theta})$.

For AI: Sufficient statistics are the foundation of exponential families, which include Gaussian, Bernoulli, Poisson, Gamma, and most common distributions. The natural parameterisation of exponential families (used in generalised linear models and variational inference) exploits sufficient statistics to simplify computation. In the VAE (Kingma & Welling, 2014), the encoder $q_\phi(\mathbf{z}|\mathbf{x})$ outputs sufficient statistics $(\boldsymbol{\mu}_\phi(\mathbf{x}), \boldsymbol{\sigma}_\phi(\mathbf{x}))$ of the Gaussian posterior approximation.

3.5 Completeness and the Lehmann-Scheffe Theorem

Definition 3.5 (Complete Statistic). A sufficient statistic $T$ is complete if for any measurable function $g$: $\mathbb{E}[g(T); \theta] = 0$ for all $\theta \in \Theta$ implies $g(T) = 0$ a.s. for all $\theta$.

Completeness means the only function of $T$ that has zero expectation everywhere is the zero function - there are no "hidden" unbiasedness conditions.

Lehmann-Scheffe Theorem. If $T$ is a complete sufficient statistic and $\tilde{\theta} = h(T)$ is unbiased for $\theta$, then $\tilde{\theta}$ is the unique minimum variance unbiased estimator (UMVUE) of $\theta$.

The UMVUE is the "best possible" unbiased estimator - it has the lowest variance among all unbiased estimators, and it is unique.

Example: For $\mathcal{N}(\mu, \sigma^2)$ with $\sigma^2$ known, $T = \bar{x}$ is complete sufficient, and $\bar{x}$ itself is unbiased for $\mu$. By Lehmann-Scheffe, $\bar{x}$ is the UMVUE of $\mu$ - there is no unbiased estimator with smaller variance. This also confirms that $\bar{x}$ achieves the CRB $\sigma^2/n$.

---

4. Fisher Information and the Cramer-Rao Bound

4.1 The Score Function

The score function is the fundamental quantity linking estimation and information theory.

Definition 4.1 (Score Function). For a single observation $x$ from $p(x;\theta)$, the score is:

$$s(\theta; x) = \frac{\partial}{\partial \theta} \log p(x; \theta)$$

For a sample of $n$ iid observations, the total score is:

$$S_n(\theta) = \sum_{i=1}^n \frac{\partial}{\partial \theta} \log p(x^{(i)};\theta)$$

The MLE solves $S_n(\hat{\theta}) = 0$ (the score equation).

Zero-mean property of the score:

$$\mathbb{E}\!\left[s(\theta; X)\right] = \int \frac{\partial \log p(x;\theta)}{\partial \theta} p(x;\theta)\, dx = \int \frac{\frac{\partial p(x;\theta)}{\partial \theta}}{p(x;\theta)} p(x;\theta)\, dx = \frac{\partial}{\partial \theta}\int p(x;\theta)\, dx = \frac{\partial}{\partial \theta}(1) = 0$$

under regularity conditions permitting interchange of differentiation and integration.

Interpretation: The score measures how sensitively the log-likelihood responds to perturbations of $\theta$. At the true parameter value, the score has zero mean - on average, the likelihood is at a stationary point. But it has positive variance, measuring how much information the data carries about $\theta$.

SCORE FUNCTION - GEOMETRIC INTUITION
========================================================================

  Log-likelihood \\ell(\\theta) = \\Sigma_i log p(x_i;\\theta) as a function of \\theta:

         \\ell(\\theta)
          |
          |         +-----+
          |       +-+     +-+
          |     +-+         +-+
          |---+-+             +-+----
          |                        \\theta
                    ^
                 \\theta* (true parameter)
                 slope = score = 0 at maximum

  Score s(\\theta; data) = d\\ell/d\\theta is the slope of the log-likelihood.
  High variance of score at \\theta* means the data sharply identifies
  \\theta* (the likelihood is "peaky"). Low variance means the data
  are relatively uninformative about \\theta.

========================================================================

4.2 Fisher Information

Definition 4.2 (Fisher Information). The Fisher information for parameter $\theta$ in model $p(x;\theta)$ is:

$$\mathcal{I}(\theta) = \mathbb{E}\!\left[\left(\frac{\partial \log p(X;\theta)}{\partial \theta}\right)^2\right] = \operatorname{Var}\!\left(\frac{\partial \log p(X;\theta)}{\partial \theta}\right)$$

(The second equality uses $\mathbb{E}[s] = 0$.)

Alternative form (under regularity conditions):

$$\mathcal{I}(\theta) = -\mathbb{E}\!\left[\frac{\partial^2 \log p(X;\theta)}{\partial \theta^2}\right]$$

Proof of equivalence: Differentiating $\mathbb{E}[s(\theta;X)] = 0$ with respect to $\theta$:

$$0 = \frac{\partial}{\partial\theta}\int \frac{\partial \log p}{\partial \theta} p\, dx = \int \left[\frac{\partial^2 \log p}{\partial \theta^2} p + \frac{\partial \log p}{\partial \theta} \cdot \frac{\partial p}{\partial \theta}\right] dx$$ $$= \mathbb{E}\!\left[\frac{\partial^2 \log p}{\partial \theta^2}\right] + \mathbb{E}\!\left[\left(\frac{\partial \log p}{\partial\theta}\right)^2\right]$$

Therefore $\mathcal{I}(\theta) = -\mathbb{E}\!\left[\frac{\partial^2 \log p}{\partial\theta^2}\right]$. $\square$

The second-derivative form has intuitive content: $\mathcal{I}(\theta)$ is the expected (negative) curvature of the log-likelihood. High curvature = the likelihood has a sharp peak = the data strongly identifies $\theta$. Low curvature = flat likelihood = data is relatively uninformative.

Fisher information for an iid sample of $n$ observations:

$$\mathcal{I}_n(\theta) = n \cdot \mathcal{I}_1(\theta)$$

Information accumulates linearly with sample size - doubling data doubles information, halving uncertainty (in variance terms: $1/\mathcal{I}_n = (1/\mathcal{I}_1)/n$).

Computed examples:

Bernoulli$(p)$: $\log p(x;p) = x\log p + (1-x)\log(1-p)$. Score: $s = x/p - (1-x)/(1-p)$. Score squared: $s^2 = \frac{(x-p)^2}{p^2(1-p)^2}$.

$$\mathcal{I}(p) = \mathbb{E}[s^2] = \frac{\operatorname{Var}(X)}{p^2(1-p)^2} = \frac{p(1-p)}{p^2(1-p)^2} = \frac{1}{p(1-p)}$$

$\mathcal{N}(\mu, \sigma^2)$ with $\sigma^2$ known: $\log p(x;\mu) = -\frac{(x-\mu)^2}{2\sigma^2} + \text{const}$. Score: $s = (x-\mu)/\sigma^2$.

$$\mathcal{I}(\mu) = \mathbb{E}\!\left[\frac{(X-\mu)^2}{\sigma^4}\right] = \frac{\operatorname{Var}(X)}{\sigma^4} = \frac{\sigma^2}{\sigma^4} = \frac{1}{\sigma^2}$$

Poisson$(\lambda)$: $\log p(x;\lambda) = x\log\lambda - \lambda - \log(x!)$. Score: $s = x/\lambda - 1$.

$$\mathcal{I}(\lambda) = \mathbb{E}\!\left[\left(\frac{X}{\lambda} - 1\right)^2\right] = \frac{\operatorname{Var}(X)}{\lambda^2} = \frac{\lambda}{\lambda^2} = \frac{1}{\lambda}$$

Pattern: For exponential family distributions, $\mathcal{I}(\theta) = 1/\operatorname{Var}(\hat{\theta}_{\text{MLE}})$ in scalar cases - information is the reciprocal of the natural parameter variance.

Multivariate Fisher Information Matrix. For $\boldsymbol{\theta} \in \mathbb{R}^k$:

$$\mathcal{I}(\boldsymbol{\theta})_{jl} = \mathbb{E}\!\left[\frac{\partial \log p}{\partial \theta_j} \cdot \frac{\partial \log p}{\partial \theta_l}\right] = -\mathbb{E}\!\left[\frac{\partial^2 \log p}{\partial \theta_j \partial \theta_l}\right]$$

In matrix form: $\mathcal{I}(\boldsymbol{\theta}) = \mathbb{E}[\mathbf{s}\mathbf{s}^\top] = -\mathbb{E}[H_{\log p}]$ where $H_{\log p}$ is the Hessian of the log-likelihood.

The FIM is always positive semi-definite (as an outer product expectation), and positive definite under identifiability.

4.3 The Cramer-Rao Lower Bound

The Cramer-Rao Lower Bound (CRB) is the fundamental limit on estimation accuracy.

Theorem 4.3 (Cramer-Rao Lower Bound). Let $\hat{\theta}$ be any unbiased estimator of $\theta$ based on $n$ iid observations. Under regularity conditions:

$$\operatorname{Var}(\hat{\theta}) \geq \frac{1}{n\mathcal{I}(\theta)}$$

For the multivariate case with unbiased estimator $\hat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta} \in \mathbb{R}^k$:

$$\operatorname{Cov}(\hat{\boldsymbol{\theta}}) \succeq \frac{1}{n}\mathcal{I}(\boldsymbol{\theta})^{-1}$$

(meaning $\operatorname{Cov}(\hat{\boldsymbol{\theta}}) - \frac{1}{n}\mathcal{I}(\boldsymbol{\theta})^{-1}$ is positive semidefinite).

Proof (scalar case, single observation):

We need to show $\operatorname{Var}(\hat{\theta}) \cdot \mathcal{I}(\theta) \geq 1$.

Since $\hat{\theta}$ is unbiased: $\mathbb{E}[\hat{\theta}] = \theta$. Differentiating with respect to $\theta$:

$$1 = \frac{\partial}{\partial\theta}\int \hat{\theta}(x)\, p(x;\theta)\, dx = \int \hat{\theta}(x)\, \frac{\partial \log p(x;\theta)}{\partial\theta}\, p(x;\theta)\, dx = \mathbb{E}[\hat{\theta} \cdot s(\theta;X)]$$

Now apply the Cauchy-Schwarz inequality to $\mathbb{E}[\hat{\theta} \cdot s]$:

$$1 = \mathbb{E}[\hat{\theta} \cdot s] = \mathbb{E}[(\hat{\theta} - \theta) \cdot s] \quad (\text{since } \mathbb{E}[s] = 0)$$

By Cauchy-Schwarz: $\mathbb{E}[(\hat{\theta} - \theta) \cdot s]^2 \leq \mathbb{E}[(\hat{\theta}-\theta)^2] \cdot \mathbb{E}[s^2]$

$$1 \leq \operatorname{Var}(\hat{\theta}) \cdot \mathcal{I}(\theta)$$ $$\operatorname{Var}(\hat{\theta}) \geq \frac{1}{\mathcal{I}(\theta)} \qquad \square$$

When is the bound tight? The CRB is achieved (with equality in Cauchy-Schwarz) if and only if $(\hat{\theta} - \theta) = c(\theta) \cdot s(\theta; X)$ for some function $c(\theta)$. This occurs precisely for exponential family distributions, and the efficient estimator is the MLE.

Biased estimator CRB. For biased estimators with $\mathbb{E}[\hat{\theta}] = \theta + b(\theta)$:

$$\operatorname{Var}(\hat{\theta}) \geq \frac{(1 + b'(\theta))^2}{\mathcal{I}(\theta)}$$

This shows that introducing bias (via regularisation) can actually reduce the CRB - a biased estimator faces a different, potentially less demanding bound.

Examples of efficient estimators:

The sample mean is efficient for all these single-parameter exponential families. This is not a coincidence - it follows from the structure of exponential families.

4.4 Fisher Information Matrix in Practice

For a neural network with parameters $\boldsymbol{\theta} \in \mathbb{R}^p$ defining a conditional distribution $p(y|\mathbf{x};\boldsymbol{\theta})$, the FIM is:

$$\mathcal{I}(\boldsymbol{\theta}) = \mathbb{E}_{(\mathbf{x},y)\sim p_{\text{data}}}\!\left[\nabla_{\boldsymbol{\theta}} \log p(y|\mathbf{x};\boldsymbol{\theta})\, \nabla_{\boldsymbol{\theta}} \log p(y|\mathbf{x};\boldsymbol{\theta})^\top\right]$$

This is an outer product of gradient vectors, making it a $p \times p$ PSD matrix. For modern LLMs with $p \sim 10^{11}$ parameters, this matrix is completely intractable to store (would require $10^{22}$ bytes).

Empirical FIM. In practice, approximate by sampling:

$$\hat{\mathcal{I}}(\boldsymbol{\theta}) = \frac{1}{n}\sum_{i=1}^n \nabla_{\boldsymbol{\theta}} \log p(y^{(i)}|\mathbf{x}^{(i)};\boldsymbol{\theta})\, \nabla_{\boldsymbol{\theta}} \log p(y^{(i)}|\mathbf{x}^{(i)};\boldsymbol{\theta})^\top$$

This is the average squared gradient - exactly what is accumulated in gradient variance estimates.

Natural gradient (Amari, 1998). The ordinary gradient $\nabla_{\boldsymbol{\theta}} \mathcal{L}$ is the steepest ascent direction in Euclidean parameter space $\mathbb{R}^p$. But parameters do not live in Euclidean space - they live in the statistical manifold of distributions, where the natural metric is the Fisher information. The natural gradient is:

$$\tilde{\nabla}_{\boldsymbol{\theta}} \mathcal{L} = \mathcal{I}(\boldsymbol{\theta})^{-1} \nabla_{\boldsymbol{\theta}} \mathcal{L}$$

This is the steepest ascent direction in the metric defined by the FIM. Natural gradient is invariant to reparametrisation: if we change from $\boldsymbol{\theta}$ to $\boldsymbol{\phi} = g(\boldsymbol{\theta})$, the natural gradient update gives the same parameter change.

K-FAC approximation (Martens & Grosse, 2015). Full FIM inversion costs $O(p^3)$. K-FAC approximates $\mathcal{I}^{-1}$ by assuming layer-wise independence and Kronecker factorisation of each layer's FIM: $\mathcal{I}^{[l]} \approx A^{[l]} \otimes G^{[l]}$, where $A^{[l]} = \mathbb{E}[\mathbf{a}^{[l-1]}(\mathbf{a}^{[l-1]})^\top]$ (input covariance) and $G^{[l]} = \mathbb{E}[\delta^{[l]}(\delta^{[l]})^\top]$ (gradient covariance). Inverting the Kronecker product: $(A \otimes G)^{-1} = A^{-1} \otimes G^{-1}$, reducing cost from $O(p^3)$ to $O(d_{\text{in}}^3 + d_{\text{out}}^3)$ per layer.

Jeffreys prior. In Bayesian statistics (preview of Section04), the Jeffreys prior $p(\theta) \propto \sqrt{\mathcal{I}(\theta)}$ is the "uninformative" prior that is invariant to reparametrisation. For $\theta = p$ in Bernoulli: $\mathcal{I}(p) = 1/(p(1-p))$, so $p_J(p) \propto p^{-1/2}(1-p)^{-1/2}$, which is $\operatorname{Beta}(1/2, 1/2)$.

Preview: MAP Estimation. Adding a prior $p(\boldsymbol{\theta})$ and finding $\arg\max_{\boldsymbol{\theta}} [p(\mathcal{D}|\boldsymbol{\theta}) p(\boldsymbol{\theta})]$ gives the Maximum A Posteriori (MAP) estimate - MLE regularised by the log-prior. MAP with a Gaussian prior $\mathcal{N}(0, \lambda^{-1}I)$ gives L2-regularised MLE. The full Bayesian treatment of priors, posteriors, and conjugate families is in Section04 Bayesian Inference.

---

5. Maximum Likelihood Estimation

5.1 The Likelihood Principle

Definition 5.1 (Likelihood Function). Given observed data $\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}$ and a parametric model $p(\cdot;\boldsymbol{\theta})$, the likelihood function is:

$$L(\boldsymbol{\theta}; \mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)}) = \prod_{i=1}^n p(\mathbf{x}^{(i)};\boldsymbol{\theta})$$

The log-likelihood is:

$$\ell(\boldsymbol{\theta}) = \log L(\boldsymbol{\theta}) = \sum_{i=1}^n \log p(\mathbf{x}^{(i)};\boldsymbol{\theta})$$

The Maximum Likelihood Estimator (MLE) is:

$$\hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\theta} \in \Theta} \ell(\boldsymbol{\theta}) = \arg\max_{\boldsymbol{\theta} \in \Theta} \sum_{i=1}^n \log p(\mathbf{x}^{(i)};\boldsymbol{\theta})$$

The likelihood is not a probability. $L(\boldsymbol{\theta}; \mathbf{x})$ is a function of $\boldsymbol{\theta}$ for fixed data $\mathbf{x}$, not a probability distribution over $\boldsymbol{\theta}$. It does not integrate to 1 over $\boldsymbol{\theta}$. The statement "$L(\boldsymbol{\theta}) = 0.3$" means "the data has probability density 0.3 under parameter $\boldsymbol{\theta}$", not "there is a 30% chance $\boldsymbol{\theta}$ is the true parameter".

Why log? Three reasons:

1. Computational: $\prod_{i=1}^n p_i$ underflows to zero for large $n$ (e.g., $0.5^{1000} \approx 10^{-301}$); sums are numerically stable
2. Mathematical: sums are easier to differentiate and optimise than products
3. Statistical: log is a monotone transformation, so $\arg\max L = \arg\max \log L$

The likelihood principle states that all information about $\boldsymbol{\theta}$ in the data is contained in the likelihood function. Two datasets with proportional likelihoods should lead to the same inferences about $\boldsymbol{\theta}$. This is a philosophical principle - frequentists may not accept it fully (confidence intervals depend on the sampling procedure, not just the likelihood), but it underlies MLE and Bayesian inference alike.

5.2 Deriving MLEs: Scalar Parameters

The standard procedure for finding the MLE:

1. Write the log-likelihood $\ell(\theta) = \sum_i \log p(x^{(i)};\theta)$
2. Differentiate and set $\ell'(\theta) = 0$ (the score equation)
3. Verify it is a maximum (second derivative $< 0$ or global structure)
4. Check boundary cases (the maximum might be at $\Theta$'s boundary)

Bernoulli MLE. Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} \operatorname{Bern}(p)$.

$$\ell(p) = \sum_{i=1}^n \left[x^{(i)}\log p + (1-x^{(i)})\log(1-p)\right] = S\log p + (n-S)\log(1-p)$$

where $S = \sum_{i=1}^n x^{(i)}$ is the total number of successes. Setting $\ell'(p) = 0$:

$$\frac{S}{p} - \frac{n-S}{1-p} = 0 \implies S(1-p) = p(n-S) \implies \hat{p}_{\text{MLE}} = \frac{S}{n} = \bar{x}$$

The MLE is the sample proportion - the intuitive answer.

Poisson MLE. Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} \operatorname{Poisson}(\lambda)$.

$$\ell(\lambda) = \sum_{i=1}^n \left[x^{(i)}\log\lambda - \lambda - \log(x^{(i)}!)\right] = S\log\lambda - n\lambda + \text{const}$$

Setting $\ell'(\lambda) = S/\lambda - n = 0$: $\hat{\lambda}_{\text{MLE}} = S/n = \bar{x}$.

The MLE of the Poisson rate is the sample mean - the estimator of the mean equals the MLE because $\mathbb{E}[X] = \lambda$.

Exponential MLE. Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} \operatorname{Exp}(\lambda)$, i.e., $p(x;\lambda) = \lambda e^{-\lambda x}$ for $x > 0$.

$$\ell(\lambda) = n\log\lambda - \lambda\sum_{i=1}^n x^{(i)}$$

Setting $\ell'(\lambda) = n/\lambda - \sum x^{(i)} = 0$: $\hat{\lambda}_{\text{MLE}} = n/\sum x^{(i)} = 1/\bar{x}$.

The MLE of the rate is the reciprocal of the sample mean - since $\mathbb{E}[X] = 1/\lambda$, this is natural.

Uniform MLE. Let $x^{(1)}, \ldots, x^{(n)} \overset{\text{iid}}{\sim} \mathcal{U}(0,\theta)$, $p(x;\theta) = 1/\theta$ for $x \in [0,\theta]$.

$$\ell(\theta) = -n\log\theta \cdot \mathbb{1}[\theta \geq \max_i x^{(i)}]$$

The likelihood is zero for $\theta < \max_i x^{(i)}$ (some observation would be outside $[0,\theta]$) and decreasing in $\theta$ for $\theta \geq \max_i x^{(i)}$. The MLE is at the boundary: $\hat{\theta}_{\text{MLE}} = \max_i x^{(i)}$.

This is an example where the score equation gives no solution (the log-likelihood has no interior critical point); the MLE is found by reasoning about the likelihood's shape. Also, the MLE here is biased: $\mathbb{E}[\max_i x^{(i)}] = n\theta/(n+1) < \theta$.

5.3 Deriving MLEs: Multivariate Gaussian

Let $\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(n)} \overset{\text{iid}}{\sim} \mathcal{N}(\boldsymbol{\mu}, \Sigma)$ with $\boldsymbol{\mu} \in \mathbb{R}^d$ and $\Sigma \in \mathbb{S}^d_{++}$.

The log-likelihood is:

$$\ell(\boldsymbol{\mu}, \Sigma) = -\frac{n}{2}\log\det(\Sigma) - \frac{1}{2}\sum_{i=1}^n (\mathbf{x}^{(i)} - \boldsymbol{\mu})^\top \Sigma^{-1}(\mathbf{x}^{(i)} - \boldsymbol{\mu}) - \frac{nd}{2}\log(2\pi)$$

MLE of $\boldsymbol{\mu}$: Taking the gradient with respect to $\boldsymbol{\mu}$ and setting to zero:

$$\nabla_{\boldsymbol{\mu}} \ell = \sum_{i=1}^n \Sigma^{-1}(\mathbf{x}^{(i)} - \boldsymbol{\mu}) = 0 \implies \hat{\boldsymbol{\mu}}_{\text{MLE}} = \frac{1}{n}\sum_{i=1}^n \mathbf{x}^{(i)} = \bar{\mathbf{x}}$$

MLE of $\Sigma$: Substituting $\hat{\boldsymbol{\mu}} = \bar{\mathbf{x}}$ and differentiating with respect to $\Sigma^{-1}$ (using matrix calculus identities $\nabla_A \log\det A = A^{-\top}$ and $\nabla_A \operatorname{tr}(BA) = B^\top$):

$$\frac{\partial \ell}{\partial \Sigma^{-1}} = \frac{n}{2}\Sigma - \frac{1}{2}\sum_{i=1}^n (\mathbf{x}^{(i)} - \bar{\mathbf{x}})(\mathbf{x}^{(i)} - \bar{\mathbf{x}})^\top = 0$$ $$\hat{\Sigma}_{\text{MLE}} = \frac{1}{n}\sum_{i=1}^n (\mathbf{x}^{(i)} - \bar{\mathbf{x}})(\mathbf{x}^{(i)} - \bar{\mathbf{x}})^\top$$

Bias of $\hat{\Sigma}_{\text{MLE}}$: By the same calculation as in the scalar case, $\mathbb{E}[\hat{\Sigma}_{\text{MLE}}] = \frac{n-1}{n}\Sigma$. The unbiased estimator is $S = \frac{1}{n-1}\sum_i (\mathbf{x}^{(i)} - \bar{\mathbf{x}})(\mathbf{x}^{(i)} - \bar{\mathbf{x}})^\top$.

For AI: Fitting a Gaussian to activations or embeddings is done in numerous ML methods. The empirical covariance of a layer's activations, used in whitening transforms, weight initialisation (Xavier/He initialisation matches the second moment), and covariance-regularised fine-tuning, all use this MLE formula. The multivariate Gaussian MLE is also the building block for Gaussian mixture models (EM algorithm) and linear discriminant analysis.

5.4 MLE as Cross-Entropy Minimisation

This connection is perhaps the most important in all of applied ML.

Theorem 5.4. For a classification model $p(y|\mathbf{x};\boldsymbol{\theta})$ and iid data $\{(\mathbf{x}^{(i)}, y^{(i)})\}_{i=1}^n$:

$$\hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\theta}} \sum_{i=1}^n \log p(y^{(i)}|\mathbf{x}^{(i)};\boldsymbol{\theta}) = \arg\min_{\boldsymbol{\theta}} \left[-\frac{1}{n}\sum_{i=1}^n \log p(y^{(i)}|\mathbf{x}^{(i)};\boldsymbol{\theta})\right]$$

The right-hand side is the negative log-likelihood (NLL) or cross-entropy loss:

$$\mathcal{L}_{\text{CE}}(\boldsymbol{\theta}) = -\frac{1}{n}\sum_{i=1}^n \log p(y^{(i)}|\mathbf{x}^{(i)};\boldsymbol{\theta})$$

Connection to KL divergence: Define the empirical distribution $p_{\text{data}}(y|\mathbf{x}) = \frac{1}{n}\sum_i \delta_{(x^{(i)}, y^{(i)})}$. Then:

$$\mathcal{L}_{\text{CE}} = -\frac{1}{n}\sum_i \log p(y^{(i)}|\mathbf{x}^{(i)};\boldsymbol{\theta}) = H(p_{\text{data}}) + D_{\mathrm{KL}}(p_{\text{data}} \| p_\theta)$$

Since $H(p_{\text{data}})$ does not depend on $\boldsymbol{\theta}$:

$$\arg\min_{\boldsymbol{\theta}} \mathcal{L}_{\text{CE}} = \arg\min_{\boldsymbol{\theta}} D_{\mathrm{KL}}(p_{\text{data}} \| p_{\boldsymbol{\theta}})$$

MLE minimises the KL divergence from the data distribution to the model. This is the information-theoretic interpretation of MLE.

For language models: An LLM defines $p_{\boldsymbol{\theta}}(x_t | x_1, \ldots, x_{t-1})$. Training by NLL is MLE:

$$\hat{\boldsymbol{\theta}} = \arg\min_{\boldsymbol{\theta}} -\frac{1}{T}\sum_{t=1}^T \log p_{\boldsymbol{\theta}}(x_t | x_{<t})$$

where the sum is over all tokens in the training corpus. This objective is used verbatim in every large language model - GPT-4 (OpenAI, 2023), LLaMA-3 (Meta, 2024), Gemini (Google, 2023), Claude (Anthropic, 2024). The perplexity $\exp(\mathcal{L}_{\text{CE}})$ is the standard evaluation metric.

Label smoothing as regularised MLE. Standard MLE uses hard labels $y^{(i)} \in \{0,1\}^K$ (one-hot). Label smoothing (Szegedy et al., 2016) uses soft targets $\tilde{y}^{(i)} = (1-\varepsilon)y^{(i)} + \varepsilon/K$, adding a small weight to incorrect classes. This is equivalent to regularising MLE by mixing the empirical distribution with a uniform prior - reducing overconfidence and improving calibration.

5.5 Properties of MLE

Under regularity conditions (the "Cramer-Rao conditions" on the model), MLE satisfies:

1. Consistency: $\hat{\boldsymbol{\theta}}_{\text{MLE}} \xrightarrow{P} \boldsymbol{\theta}^*$. The proof uses the fact that the expected log-likelihood $\mathbb{E}[\ell(\boldsymbol{\theta})/n]$ is uniquely maximised at the true parameter $\boldsymbol{\theta}^*$ (by the strict convexity of KL divergence), and the LLN ensures $\ell(\boldsymbol{\theta})/n \to \mathbb{E}[\log p(X;\boldsymbol{\theta})]$ uniformly.

2. Asymptotic normality:

$$\sqrt{n}\left(\hat{\boldsymbol{\theta}}_{\text{MLE}} - \boldsymbol{\theta}^*\right) \xrightarrow{d} \mathcal{N}\!\left(\mathbf{0},\, \mathcal{I}(\boldsymbol{\theta}^*)^{-1}\right)$$

For large $n$: $\hat{\boldsymbol{\theta}}_{\text{MLE}} \approx \mathcal{N}(\boldsymbol{\theta}^*, \frac{1}{n}\mathcal{I}(\boldsymbol{\theta}^*)^{-1})$.

3. Asymptotic efficiency: The asymptotic covariance $\frac{1}{n}\mathcal{I}^{-1}$ achieves the CRB. No consistent estimator can have smaller asymptotic variance (at every $\boldsymbol{\theta}$). MLE is asymptotically optimal.

4. Invariance under reparametrisation: If $\hat{\boldsymbol{\theta}}$ is the MLE of $\boldsymbol{\theta}$, then $g(\hat{\boldsymbol{\theta}})$ is the MLE of $g(\boldsymbol{\theta})$ for any function $g$ (not required to be injective). This is the invariance property of MLE.

Examples of invariance:

• MLE of $\sigma$ is $\hat{\sigma} = \sqrt{\hat{\sigma}^2_{\text{MLE}}}$ (not $s$, which is the unbiased estimator)
• MLE of $1/\lambda$ in Exponential is $\bar{x}$ (consistent with $\hat{\lambda} = 1/\bar{x}$)
• MLE of $p^2$ in Bernoulli is $\hat{p}^2 = \bar{x}^2$ (biased, but the MLE)

Regularity conditions (Cramer-Rao conditions). These are the sufficient conditions for the above properties:

1. The parameter space $\Theta$ is an open subset of $\mathbb{R}^k$
2. The model is identifiable
3. The support $\{x: p(x;\boldsymbol{\theta}) > 0\}$ does not depend on $\boldsymbol{\theta}$ (excludes Uniform)
4. The log-likelihood is three times differentiable in $\boldsymbol{\theta}$
5. The Fisher information matrix is positive definite at $\boldsymbol{\theta}^*$

When these fail (as in the Uniform example), MLE may still be the natural estimator but its properties differ.

5.6 Numerical MLE

For complex models, the score equation $\ell'(\boldsymbol{\theta}) = 0$ has no closed-form solution and must be solved numerically.

Gradient ascent on log-likelihood. The simplest approach:

$$\boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t + \eta \nabla_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta}_t)$$

For neural networks trained with mini-batches: stochastic gradient ascent on the mini-batch log-likelihood. This is exactly SGD on NLL loss with a negated gradient.

Newton-Raphson / Fisher scoring. Use second-order information:

$$\boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t - [H_\ell(\boldsymbol{\theta}_t)]^{-1} \nabla \ell(\boldsymbol{\theta}_t)$$

where $H_\ell$ is the Hessian of the log-likelihood. Replacing the Hessian with the expected negative Hessian (Fisher information) gives Fisher scoring:

$$\boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t + \frac{1}{n}\mathcal{I}(\boldsymbol{\theta}_t)^{-1} \nabla \ell(\boldsymbol{\theta}_t)$$

This is exactly the natural gradient update. Newton-Raphson converges quadratically near the maximum - much faster than gradient ascent - but requires inverting a $k \times k$ matrix per step.

EM algorithm. For models with latent variables (mixtures, HMMs, VAEs), the EM algorithm alternates between:

• E-step: compute $Q(\boldsymbol{\theta};\boldsymbol{\theta}_t) = \mathbb{E}_{z|\mathbf{x};\boldsymbol{\theta}_t}[\log p(\mathbf{x},z;\boldsymbol{\theta})]$
• M-step: $\boldsymbol{\theta}_{t+1} = \arg\max_{\boldsymbol{\theta}} Q(\boldsymbol{\theta};\boldsymbol{\theta}_t)$

EM guarantees that $\ell(\boldsymbol{\theta}_{t+1}) \geq \ell(\boldsymbol{\theta}_t)$ - the log-likelihood is non-decreasing. EM is a specialised case of variational inference (to be developed in Section04 Bayesian Inference).

Numerical pitfalls:

• Overflow/underflow: compute $\ell = \sum \log p_i$, not $\log \prod p_i$ (avoid multiplying probabilities)
• Log-sum-exp trick: for softmax-based likelihoods, $\log \sum_k e^{z_k} = z_{\max} + \log\sum_k e^{z_k - z_{\max}}$
• Flat log-likelihoods: when $\mathcal{I}(\boldsymbol{\theta}) \approx 0$ near the solution, gradient methods converge extremely slowly; natural gradient methods help
• Local maxima: the log-likelihood may be multimodal for mixture models and neural networks; multiple restarts are necessary

---