%matplotlib inline
from d2l import torch as d2l
import numpy as onp
import torch26.3 Maximum Likelihood
Many standard machine-learning losses have a probabilistic interpretation. When we minimize a classifier’s cross-entropy or a regressor’s squared error under a Gaussian noise model, we are answering one question: which parameters make the observed data most probable? That is the maximum likelihood principle, the idea that turns a probabilistic model into a trainable objective. This section states the principle and proves the three equivalences that make it operational: likelihood becomes a negative log-likelihood, the negative log-likelihood is the cross-entropy to the data, and the Gaussian case is mean squared error. It then asks how good the resulting estimates are (consistency, asymptotic normality, and the Fisher information) and shows how putting a prior back in recovers \(L_2\) regularization. The section ends where modern generative models begin: when latent variables make the likelihood intractable, a lower bound (the ELBO) and the EM algorithm make it tractable again. We keep the running coin-flip example throughout because it is small enough to solve by hand, which lets us check every claim against the answer we already know.
%matplotlib inline
from d2l import tensorflow as d2l
import numpy as onp
import tensorflow as tf%matplotlib inline
from d2l import jax as d2l
import jax
from jax import numpy as jnp
import numpy as onp%matplotlib inline
from d2l import mxnet as d2l
from mxnet import autograd, np, npx
import numpy as onp
npx.set_np()26.3.1 The Maximum Likelihood Principle
Suppose we have a model with parameters \(\boldsymbol{\theta}\) and observed data \(X\). The likelihood is the model probability or density of the observed data, viewed as a function of the parameter:
\[ L(\boldsymbol{\theta};X)=P(X\mid\boldsymbol{\theta}). \]
Maximum likelihood is the direct rule
\[ \hat{\boldsymbol{\theta}}_{\textrm{MLE}} = \mathop{\mathrm{argmax}}_{\boldsymbol{\theta}} P(X \mid \boldsymbol{\theta}). \tag{26.3.1}\]
No prior is part of this definition. MLE asks which member of the chosen model family fits the observations best; it does not assign probabilities to the parameters themselves. Later, Section 26.3.4 introduces the distinct maximum a posteriori rule \(\operatorname*{argmax}_{\boldsymbol{\theta}}P(\boldsymbol{\theta}\mid X)\), which multiplies the likelihood by a prior density. A prior that is constant in the chosen parameterization makes the two optimizers coincide, but “flat” is not invariant under nonlinear reparameterization and is not the foundation of maximum likelihood.
26.3.1.1 A Worked Example: The Coin
Take a single parameter \(\theta\), the probability of heads, so tails has probability \(1-\theta\). Because the flips are independent, their probabilities multiply: if \(X\) has \(n_H\) heads and \(n_T\) tails,
\[ P(X \mid \theta) = \theta^{n_H}(1-\theta)^{n_T}. \tag{26.3.2}\]
Flip \(13\) coins and observe the sequence “HHHTHTTHHHHHT”, which has \(n_H=9\) and \(n_T=4\), giving the likelihood \(P(X \mid \theta) = \theta^9(1-\theta)^4\). We know the answer we want: asked “\(9\) of \(13\) flips were heads, what is the bias?”, everyone replies \(9/13\). The value of the maximum-likelihood machinery is that it derives that number from first principles, in a way that scales to models with billions of parameters where no intuition is available.
The likelihood is a function of \(\theta\) alone, so we can plot it and read off its peak.
import numpy as onp
theta = onp.arange(0, 1, 0.001)
p = theta**9 * (1 - theta)**4.
d2l.plot(theta, p, 'theta', 'likelihood')
The curve peaks near the expected \(9/13 \approx 0.69\). To pin the maximum exactly we use calculus: at an interior maximum the derivative vanishes. Setting \(\frac{d}{d\theta}P(X\mid\theta)=0\),
\[ 0 = \frac{d}{d\theta}\,\theta^9(1-\theta)^4 = 9\theta^8(1-\theta)^4 - 4\theta^9(1-\theta)^3 = \theta^8(1-\theta)^3(9-13\theta). \]
The roots \(\theta=0\) and \(\theta=1\) assign probability zero to a sequence that contains both heads and tails, so they are minima. The remaining root is the maximum likelihood estimate,
\[ \hat\theta = \frac{9}{13}, \]
exactly the observed fraction of heads. The same calculation on the general likelihood Equation 26.3.2 gives \(\hat\theta = n_H/(n_H+n_T)\): the MLE of a coin’s bias is always the empirical frequency of heads, our first sign that maximum likelihood recovers the “obvious” estimator.
26.3.1.2 The Negative Log-Likelihood
The closed-form trick above does not extend to a real model: with billions of parameters there is no polynomial to factor, and we optimize numerically instead. But the likelihood itself is unusable in floating point. With \(n\) independent data points it is a product of \(n\) probabilities,
\[ P(X\mid\boldsymbol{\theta}) = \prod_{i=1}^n p(x_i\mid\boldsymbol{\theta}), \]
each in \([0,1]\). A billion factors near \(1/2\) produce roughly \((1/2)^{10^9}\), underflowing to zero in any floating-point format long before we can differentiate it.
The logarithm rescues us, turning the product into a sum that stays in range:
\[ \log\!\big((1/2)^{10^9}\big) = 10^9 \cdot \log(1/2) \approx -6.93\times10^{8}, \]
which fits comfortably even in single precision. Since \(x \mapsto \log x\) is strictly increasing, it preserves the location of the maximum, so maximizing the likelihood is the same as maximizing the log-likelihood \(\log P(X\mid\boldsymbol{\theta})\). We use this exact reasoning for the naive Bayes classifier in Section 26.7. Finally, because we prefer to minimize losses, we flip the sign and define the negative log-likelihood (NLL):
\[ \ell(\boldsymbol{\theta}) = -\log P(X\mid\boldsymbol{\theta}) = -\sum_{i=1}^n \log p(x_i\mid\boldsymbol{\theta}). \tag{26.3.3}\]
For the coin this reads \(\ell(\theta) = -\big(n_H\log\theta + n_T\log(1-\theta)\big)\), which we can minimize by gradient descent even for billions of flips, no closed form required. The cell below does exactly this for almost nine million flips and confirms it recovers \(n_H/(n_H+n_T)\).
# Set up our data
n_H = 8675309
n_T = 256245
# Initialize our parameters
theta = torch.tensor(0.5, requires_grad=True)
# Perform gradient descent
lr = 1e-9
for iter in range(100):
loss = -(n_H * torch.log(theta) + n_T * torch.log(1 - theta))
loss.backward()
with torch.no_grad():
theta -= lr * theta.grad
theta.grad.zero_()
# Check output
theta, n_H / (n_H + n_T)(tensor(0.9713, requires_grad=True), 0.9713101437890875)# Set up our data
n_H = 8675309
n_T = 256245
# Initialize our parameters
theta = tf.Variable(tf.constant(0.5))
# Perform gradient descent
lr = 1e-9
for iter in range(100):
with tf.GradientTape() as t:
loss = -(n_H * tf.math.log(theta) + n_T * tf.math.log(1 - theta))
theta.assign_sub(lr * t.gradient(loss, theta))
# Check output
theta, n_H / (n_H + n_T)(<tf.Variable 'Variable:0' shape=() dtype=float32, numpy=0.9713100790977478>,
0.9713101437890875)# Set up our data
n_H = 8675309
n_T = 256245
# Initialize our parameters
theta = jnp.float32(0.5)
# Define loss function
def nll(theta):
return -(n_H * jnp.log(theta) + n_T * jnp.log(1 - theta))
grad_fn = jax.grad(nll)
# Perform gradient descent
lr = 1e-9
for iter in range(100):
theta = theta - lr * grad_fn(theta)
# Check output
theta, n_H / (n_H + n_T)(Array(0.9713101, dtype=float32), 0.9713101437890875)# Set up our data
n_H = 8675309
n_T = 256245
# Initialize our parameters
theta = np.array(0.5)
theta.attach_grad()
# Perform gradient descent
lr = 1e-9
for iter in range(100):
with autograd.record():
loss = -(n_H * np.log(theta) + n_T * np.log(1 - theta))
loss.backward()
theta -= lr * theta.grad
# Check output
theta, n_H / (n_H + n_T)[14:02:32] /home/smola/mxnet/src/base.cc:48: GPU context requested, but no GPUs found.(array(0.9713101), 0.9713101437890875)The learning rate \(10^{-9}\) looks absurd. It works only because this loss is a sum over nearly nine million flips, so its gradient carries that factor of \(n\). Shrink \(n_H\) and \(n_T\) a thousandfold and the same learning rate barely moves \(\theta\); scale them up and the update diverges. Averaging the loss over examples instead of summing (as the minibatch losses in this book do) removes the coupling and returns sensible learning rates to sensible magnitudes. The lesson generalizes: a learning rate is always tuned to a loss scale, and silently changing one without the other is a classic way to break training.
Beyond numerical stability, the log gives gradients an additive form. The product rule applied to \(\prod_i p(x_i\mid\boldsymbol{\theta})\) produces \(n\) terms, each a product of \(n-1\) surviving factors; the logarithm collapses this to a single sum,
\[ \nabla_{\boldsymbol{\theta}}\, \ell(\boldsymbol{\theta}) = -\sum_{i=1}^n \frac{\nabla_{\boldsymbol{\theta}}\, p(x_i\mid\boldsymbol{\theta})}{p(x_i\mid\boldsymbol{\theta})}, \]
a chain-rule application term by term with no shared subproducts to track. This is why every minibatch loss in the book is an average of per-example NLLs: the gradient of a sum is the sum of gradients, so we can estimate \(\nabla\ell\) from a random subset of the data.
26.3.2 Maximum Likelihood Is Minimizing a Loss
We now make precise the claim from the introduction. “Minimize the loss” and “do maximum likelihood” are the same instruction in different words. The bridge is one rescaling.
26.3.2.1 NLL Is the Cross-Entropy to the Data
Dividing the NLL Equation 26.3.3 by the number of examples turns it into an average, the natural per-example loss, and that average has an information-theoretic name. Let \(\hat p_{\textrm{data}}\) be the empirical distribution, which places mass \(1/n\) on each observation; masses at repeated values stack, so an outcome \(x\) that occurs \(n_x\) times in the sample carries total mass \(\hat p_{\textrm{data}}(x) = n_x/n\), its empirical frequency. (We work in natural logarithms, so the units are nats; switching to bits merely rescales every quantity by \(\ln 2\) and changes no \(\mathop{\mathrm{argmin}}\).)
Proposition (discrete MLE = minimum cross-entropy). For discrete outcomes, maximizing the likelihood is the same as minimizing the cross-entropy from the empirical distribution \(\hat p_{\textrm{data}}\) to the model \(p_{\boldsymbol{\theta}}\):
\[ \mathop{\mathrm{argmax}}_{\boldsymbol{\theta}} \prod_{i=1}^n p(x_i\mid\boldsymbol{\theta}) = \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}}\, \textrm{CE}\!\left(\hat p_{\textrm{data}},\, p_{\boldsymbol{\theta}}\right), \qquad \textrm{CE}(\hat p_{\textrm{data}}, p_{\boldsymbol{\theta}}) = -\!\!\sum_{x} \hat p_{\textrm{data}}(x)\log p_{\boldsymbol{\theta}}(x). \]
Proof. Take the average negative log-likelihood and regroup the sum over data points into a sum over distinct outcomes, pooling the observations with \(x_i = x\). Each observation contributes mass \(1/n\), so outcome \(x\) receives total mass \(n_x/n = \hat p_{\textrm{data}}(x)\). On the coin this regrouping is concrete: the nine head-flips pool into nine copies of \(\log p_{\boldsymbol{\theta}}(\textrm{H})\), so heads enters the average with weight \(\hat p_{\textrm{data}}(\textrm{H}) = 9/13\), and tails with weight \(4/13\). In general,
\[ \frac{1}{n}\,\ell(\boldsymbol{\theta}) = -\frac{1}{n}\sum_{i=1}^n \log p_{\boldsymbol{\theta}}(x_i) = -\sum_{x} \frac{n_x}{n}\,\log p_{\boldsymbol{\theta}}(x) = -\sum_{x} \hat p_{\textrm{data}}(x)\,\log p_{\boldsymbol{\theta}}(x) = \textrm{CE}\!\left(\hat p_{\textrm{data}},\, p_{\boldsymbol{\theta}}\right), \]
which is exactly the cross-entropy named in the proposition. Scaling by the positive constant \(1/n\) and flipping the sign do not move the optimizer, so the \(\mathop{\mathrm{argmax}}\) of the likelihood equals the \(\mathop{\mathrm{argmin}}\) of the cross-entropy. \(\blacksquare\)
This also explains why minimizing cross-entropy is the right thing to do. Cross-entropy decomposes as \(\textrm{CE}(P,Q)=H(P)+D_{\textrm{KL}}(P\|Q)\), the entropy of \(P\) plus the Kullback–Leibler divergence from \(P\) to \(Q\) (see Section 27.1 for that background). Setting \(P=\hat p_{\textrm{data}}\) and \(Q=p_{\boldsymbol{\theta}}\), the empirical entropy \(H(\hat p_{\textrm{data}})\) is fixed by the data and does not depend on \(\boldsymbol{\theta}\), so minimizing the cross-entropy is identical to minimizing the KL term alone:
\[ \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}}\, \textrm{CE}\!\left(\hat p_{\textrm{data}},\, p_{\boldsymbol{\theta}}\right) = \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}}\, D_{\textrm{KL}}\!\left(\hat p_{\textrm{data}}\,\|\,p_{\boldsymbol{\theta}}\right). \]
In other words, maximum likelihood is the KL-projection of the empirical distribution onto the model family: it picks the \(p_{\boldsymbol{\theta}}\) that is KL-closest to the data. The projection reading even survives a misspecified model: when no \(p_{\boldsymbol{\theta}}\) equals the data-generating distribution there is no true parameter to converge to, and the MLE instead converges to the parameter whose \(p_{\boldsymbol{\theta}}\) is KL-closest to the data-generating distribution (White 1982). For this fixed discrete empirical distribution, the cross-entropy cannot fall below the constant \(H(\hat p_{\textrm{data}})\); optimization changes only the KL gap, as Figure 26.3.1 shows. This is an algebraic decomposition of the observed frequencies, not by itself a statement about irreducible population randomness. As a special case, the categorical NLL of a classifier with one-hot labels is precisely the softmax cross-entropy loss of Section 3.1.
26.3.2.2 From Probabilities to Densities
Everything so far was phrased for discrete outcomes, where \(P(X\mid\boldsymbol{\theta})\) is a genuine probability. For continuous data we replace probabilities by densities \(p\). The discrete empirical-cross-entropy identity above should not be copied literally: an atomic empirical measure is generally singular with respect to a continuous density, so that KL divergence is infinite. The sample objective remains the average negative log-density, while its population expectation is a cross-entropy and differs from \(D_{\textrm{KL}}(p_\star\|p_{\boldsymbol{\theta}})\) by the fixed differential entropy of the data-generating density when those quantities are finite. The remaining practical issue is that the probability of observing any exact real value is zero, so the naive probability likelihood is \(0\) for every \(\boldsymbol{\theta}\). The resolution is to ask for a match only to within a small tolerance \(\epsilon\), then watch the \(\epsilon\) cancel.
For i.i.d. observations \(x_1,\ldots,x_n\), the probability that each lands in a window of width \(\epsilon\) is, to first order,
\[ P\big(X_i \in [x_i, x_i+\epsilon]\ \forall i \mid \boldsymbol{\theta}\big) \approx \epsilon^n \prod_{i=1}^n p(x_i\mid\boldsymbol{\theta}). \]
Taking the negative log,
\[ -\log P\big(X_i \in [x_i, x_i+\epsilon]\ \forall i\mid\boldsymbol{\theta}\big) \approx -n\log\epsilon - \sum_{i} \log p(x_i\mid\boldsymbol{\theta}). \]
The tolerance enters only through the additive constant \(-n\log\epsilon\), which is free of \(\boldsymbol{\theta}\). Demanding four digits of precision or four hundred changes that constant but never the minimizer, so we drop it and minimize
\[ -\sum_{i} \log p(x_i\mid\boldsymbol{\theta}) \]
exactly as in the discrete case. Maximum likelihood thus operates on continuous variables by swapping probabilities for densities and nothing more. The most important density to swap in is the Gaussian, which we do next.
26.3.2.3 Gaussian NLL Is Mean Squared Error
The most-used regression loss falls out of the same principle by choosing a Gaussian noise model. Suppose each target is the model’s prediction plus Gaussian noise of fixed variance, \(y_i \sim \mathcal{N}(\hat y_i, \sigma^2)\) with \(\hat y_i = f_{\boldsymbol{\theta}}(\mathbf{x}_i)\).
Proposition (Gaussian NLL = MSE). For fixed-variance Gaussian targets, the negative log-likelihood equals the mean squared error up to a constant independent of \(\boldsymbol{\theta}\):
\[ -\log \prod_{i=1}^n \mathcal{N}(y_i;\,\hat y_i,\,\sigma^2) = \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-\hat y_i)^2 + \frac{n}{2}\log(2\pi\sigma^2). \tag{26.3.4}\]
Hence \(\mathop{\mathrm{argmin}}_{\boldsymbol{\theta}}\) of this NLL is the \(\mathop{\mathrm{argmin}}\) of \(\sum_i (y_i-\hat y_i)^2\).
Proof. The Gaussian density Equation 26.2.10 is \(\mathcal{N}(y;\hat y,\sigma^2) = (2\pi\sigma^2)^{-1/2}\exp\!\big(-(y-\hat y)^2/(2\sigma^2)\big)\). Take its negative log,
\[ -\log \mathcal{N}(y_i;\hat y_i,\sigma^2) = \frac{(y_i-\hat y_i)^2}{2\sigma^2} + \frac{1}{2}\log(2\pi\sigma^2), \]
then sum over the \(n\) independent examples, which is Equation 26.3.4. The second term is the same for every \(\boldsymbol{\theta}\) (the variance is fixed), so it cannot affect the minimizer, and the leading factor \(1/(2\sigma^2)\) is a positive constant. Dropping both leaves \(\sum_i (y_i-\hat y_i)^2\), the sum of squared errors. \(\blacksquare\)
So squared-error regression is maximum likelihood under the assumption that the residuals are Gaussian, the assumption the central limit theorem makes plausible whenever the noise is a sum of many small independent effects (Section 26.2). The same recipe turns every common loss into the NLL of a chosen conditional \(p(y\mid\mathbf x)\): pick the noise model and its negative log is the loss.
| Noise model \(p(y\mid\mathbf x)\) | Negative log-likelihood \(-\log p(y\mid\mathbf x)\) | Deep-learning loss |
|---|---|---|
| Gaussian, fixed variance | \(\tfrac{1}{2\sigma^2}(y-\hat y)^2 + \textrm{const}\) | mean squared error (MSE) |
| Bernoulli | \(-\,y\log\hat p - (1-y)\log(1-\hat p)\) | binary cross-entropy (BCE) |
| Categorical (one-hot \(y\)) | \(-\sum_k y_k \log \hat p_k\) | softmax cross-entropy |
| Laplace, fixed scale | \(\tfrac{1}{b}\lvert y-\hat y\rvert + \textrm{const}\) | mean absolute error (MAE) |
Each row is one line of algebra (take the density, drop the \(\boldsymbol{\theta}\)-free constants) of exactly the kind we just did for the Gaussian. Picking a loss is picking a noise model.
26.3.3 Estimator Theory: Why Maximum Likelihood Works
How good is \(\hat{\boldsymbol{\theta}}\)? The first reassurance is the same KL picture. Drawing genuinely i.i.d. data from a true distribution \(p_{\boldsymbol{\theta}^\star}\), the average NLL is, by the law of large numbers (stated and proved in Section 26.5), an estimate of the expected cross-entropy \(\textrm{CE}(p_{\boldsymbol{\theta}^\star}, p_{\boldsymbol{\theta}})\), and the decomposition \(\textrm{CE} = H + D_{\textrm{KL}}\) shows that this population objective is minimized exactly at the truth \(\boldsymbol{\theta}=\boldsymbol{\theta}^\star\), where the KL term vanishes. “Exactly at” needs one hypothesis we have so far left implicit: the model must be identifiable, meaning distinct parameters give distinct distributions, \(\boldsymbol{\theta}\neq\boldsymbol{\theta}'\Rightarrow p_{\boldsymbol{\theta}}\neq p_{\boldsymbol{\theta}'}\); otherwise several parameters tie for the minimum and “the” true parameter is not even well defined. The hypothesis bites in deep learning: neural networks fail identifiability spectacularly, since permuting the hidden units of a layer (together with their weights) changes the parameter vector but not the function it computes, so for such models consistency can only ever be a statement about the fitted distribution \(p_{\hat{\boldsymbol{\theta}}}\), never about the parameter vector itself.
Minimizing the empirical average therefore targets the right minimizer. One caveat: the law of large numbers makes the empirical objective converge to the population one at each fixed \(\boldsymbol{\theta}\), and pointwise convergence of objectives does not by itself force their minimizers to converge. For that, the convergence must be uniform over the parameter space: the supremum over \(\boldsymbol{\theta}\) of the gap between the empirical and population objectives must go to zero. With identifiability and uniform convergence in hand, the conclusion is a genuine theorem (Wasserman 2013): for a well-specified model the MLE is consistent, converging in probability to the true parameter, \(\hat{\boldsymbol{\theta}}\xrightarrow{P}\boldsymbol{\theta}^\star\), so the probability of an error larger than any fixed tolerance goes to zero. The coin made this concrete: \(\hat\theta = n_H/(n_H+n_T) \to \theta^\star\) directly by the law of large numbers, and the Bernoulli family is identifiable, since different biases give different distributions. Consistency, asymptotic unbiasedness, and efficiency are the estimator-quality notions that Section 26.5 develops systematically; here they appear specialized to the MLE.
Consistency says only where \(\hat{\boldsymbol{\theta}}\) lands; the sharper statement is how fast and how tightly it concentrates. Under mild regularity conditions the MLE is asymptotically normal: the rescaled error converges in distribution to a Gaussian,
\[ \sqrt{n}\,\big(\hat{\boldsymbol{\theta}} - \boldsymbol{\theta}^\star\big) \;\xrightarrow{d}\; \mathcal{N}\!\big(\mathbf{0},\; I(\boldsymbol{\theta}^\star)^{-1}\big), \tag{26.3.5}\]
where \(I(\boldsymbol{\theta})\) is the Fisher information defined just below (Bishop 2006; Wasserman 2013). Here and below, “regularity conditions” means one standard package of hypotheses that we state once and take on faith: the density is smooth enough in \(\boldsymbol{\theta}\) that differentiation and integration can be interchanged (the tool of Section 24.4), the support of \(p(x\mid\boldsymbol{\theta})\) does not depend on \(\boldsymbol{\theta}\), and the true parameter lies in the interior of the parameter space (Wasserman 2013). The display says two things. First, the error shrinks at the \(1/\sqrt{n}\) rate, so the standard error of each component falls like \(1/\sqrt n\): halving it costs four times the data. Second, the limiting variance is exactly \(I(\boldsymbol{\theta}^\star)^{-1}\), and the Cramér–Rao bound stated below (Equation 26.3.7) says no unbiased estimator can have smaller variance than \(I(\boldsymbol{\theta}^\star)^{-1}/n\). The MLE attains the Cramér–Rao floor in the limit: it is asymptotically efficient, asymptotically the best unbiased estimator. The qualifiers matter. At finite \(n\) the MLE is generally only asymptotically unbiased, not exactly unbiased: the maximum-likelihood estimate of a Gaussian’s variance divides the sum of squares by \(n\) rather than \(n-1\) and so systematically underestimates \(\sigma^2\). Section 26.5 works out exactly this bias, identifying the \(n\)-divided estimator as the Gaussian MLE, when it derives the \(n-1\) correction. Maximum likelihood is optimal in the limit, not a guarantee of unbiasedness at every sample size.
Watching the theorem happen. Asymptotic normality is the section’s central limit theorem, and like the CLT proper it can be watched. The cell below returns to the coin: it draws \(20{,}000\) independent datasets of \(n=400\) flips with true bias \(\theta^\star=0.7\), computes the MLE \(\hat\theta=n_H/n\) on each, and histograms the rescaled errors \(\sqrt n\,(\hat\theta-\theta^\star)\) against the Gaussian that Equation 26.3.5 names. One scope note: the coin’s MLE is itself a sample mean, so for the coin the theorem reduces to the CLT proper; the demonstration illustrates the general statement on the one model where we can also verify it by hand. For the coin the Fisher information works out below to \(I(\theta)=1/(\theta(1-\theta))\), so the predicted limit is \(\mathcal N\bigl(0,\ \theta^\star(1-\theta^\star)\bigr)\) with variance \(0.21\), and the histogram lands on that particular Gaussian, center, width, and height. (One numerical nicety: \(\hat\theta\) lives on the lattice \(k/n\), so we offset the bin edges by half a lattice step to keep the discrete values off the edges.)
rng = onp.random.default_rng(0)
theta_star, n, reps = 0.7, 400, 20000
flips = rng.random((reps, n)) < theta_star # 20,000 datasets at once
z = onp.sqrt(n) * (flips.mean(axis=1) - theta_star) # sqrt(n)(theta_hat - theta*)
sigma2 = theta_star * (1 - theta_star) # = 1 / I(theta*)
hist, edges = onp.histogram(z, bins=40, range=(-2.025, 1.975), density=True)
x = onp.arange(-2, 2, 0.01)
d2l.plot(x, [onp.interp(x, (edges[:-1] + edges[1:]) / 2, hist),
onp.exp(-x**2 / (2 * sigma2)) / onp.sqrt(2 * onp.pi * sigma2)],
xlabel='sqrt(n) * (theta_hat - theta*)', ylabel='density',
legend=['20,000 replications', 'N(0, theta*(1-theta*))'])
26.3.3.1 Fisher Information and the Score
To make Equation 26.3.5 quantitative we need the object \(I(\boldsymbol{\theta})\) it contains. Define the score as the gradient of the log-likelihood with respect to the parameters, \(s(x;\boldsymbol{\theta}) = \nabla_{\boldsymbol{\theta}} \log p(x\mid\boldsymbol{\theta})\). The score is the per-example version of the NLL gradient we already met: setting \(\sum_i s(x_i;\boldsymbol{\theta}) = \mathbf 0\) is exactly the first-order condition that locates \(\hat{\boldsymbol{\theta}}\). At the true parameter the score has mean zero: differentiating the normalization \(\int p(x\mid\boldsymbol{\theta})\,dx = 1\) under the integral sign gives
\[ \mathbf 0 = \nabla_{\boldsymbol{\theta}} \int p(x\mid\boldsymbol{\theta})\,dx = \int \nabla_{\boldsymbol{\theta}}\, p(x\mid\boldsymbol{\theta})\,dx = \int s(x;\boldsymbol{\theta})\, p(x\mid\boldsymbol{\theta})\,dx = \mathbb{E}_{\boldsymbol{\theta}}\!\big[\,s(x;\boldsymbol{\theta})\,\big], \]
exactly the computation of Section 24.4, which also grants the interchange of derivative and integral (the first of the regularity conditions named above); the third equality is the identity \(\nabla_{\boldsymbol{\theta}}\, p = p\,\nabla_{\boldsymbol{\theta}} \log p = p\,s\). Because the score has mean zero, its variance equals its second moment, and that variance is the Fisher information (Fisher 1925),
\[ I(\boldsymbol{\theta}) = \operatorname{Var}_{\boldsymbol{\theta}}\!\big[\,s(x;\boldsymbol{\theta})\,\big] = \mathbb{E}_{\boldsymbol{\theta}}\!\big[\,s(x;\boldsymbol{\theta})\,s(x;\boldsymbol{\theta})^\top\,\big] = -\,\mathbb{E}_{\boldsymbol{\theta}}\!\big[\nabla^2_{\boldsymbol{\theta}} \log p(x\mid\boldsymbol{\theta})\big]. \tag{26.3.6}\]
The middle equality is the mean-zero property at work: \(\operatorname{Var}[s]=\mathbb E[ss^\top]\) once \(\mathbb E[s]=\mathbf 0\). The last equality is the information identity, and it follows by differentiating \(\mathbb{E}_{\boldsymbol{\theta}}[s]=\mathbf 0\) once more, with the same interchange granted. The product rule splits the derivative into two terms,
\[ \mathbf 0 = \nabla_{\boldsymbol{\theta}} \int s(x;\boldsymbol{\theta})\,p(x\mid\boldsymbol{\theta})\,dx = \int \nabla^2_{\boldsymbol{\theta}} \log p(x\mid\boldsymbol{\theta})\; p(x\mid\boldsymbol{\theta})\,dx + \int s(x;\boldsymbol{\theta})\,\big(\nabla_{\boldsymbol{\theta}}\, p(x\mid\boldsymbol{\theta})\big)^\top dx, \]
and substituting \(\nabla_{\boldsymbol{\theta}}\, p = p\,s\) in the second term turns it into \(\mathbb{E}_{\boldsymbol{\theta}}[\,s\,s^\top]\), so
\[ \mathbb{E}_{\boldsymbol{\theta}}\!\big[\nabla^2_{\boldsymbol{\theta}} \log p(x\mid\boldsymbol{\theta})\big] + \mathbb{E}_{\boldsymbol{\theta}}\!\big[\,s\,s^\top\,\big] = \mathbf 0, \]
which is the last equality of Equation 26.3.6 rearranged. The identity says that the information is the expected curvature of the NLL. The reading is intuitive: a sharply peaked log-likelihood (large curvature, large \(I\)) pins the parameter down tightly, so the estimate has small variance; a flat one leaves \(\boldsymbol{\theta}\) poorly determined. That is precisely what Equation 26.3.5 says, \(\operatorname{Var}(\hat{\boldsymbol{\theta}}) \approx I(\boldsymbol{\theta}^\star)^{-1}/n\). Figure 26.3.2 draws the picture for a Gaussian’s two parameters: the NLL is a bowl around the MLE, and the inverse of its curvature is the ellipse within which the estimate scatters.
The curvature intuition extends to a theorem constraining every estimator. The Cramér–Rao bound (Cramér 1946; Radhakrishna Rao 1945) states that any unbiased estimator \(\tilde\theta\) built from \(n\) i.i.d. observations obeys
\[ \operatorname{Var}\big(\tilde\theta\big) \;\ge\; \frac{1}{n\,I(\theta)} \tag{26.3.7}\]
(stated here for a scalar parameter; for vectors the statement reads \(\operatorname{Cov}(\tilde{\boldsymbol{\theta}}) \succeq I(\boldsymbol{\theta})^{-1}/n\) in the positive-semidefinite order) (Wasserman 2013). Information caps precision: no cleverness in constructing the estimator can beat the curvature the model itself supplies, and comparing with Equation 26.3.5 shows the MLE’s limiting variance meets the Cramér–Rao floor exactly, which is precisely the asymptotic efficiency claimed above.
Back to the coin. A single flip has log-likelihood \(\log p(x\mid\theta) = x\log\theta + (1-x)\log(1-\theta)\) for \(x\in\{0,1\}\), so the score and its derivative are
\[ s(x;\theta) = \frac{x}{\theta} - \frac{1-x}{1-\theta}, \qquad \frac{d}{d\theta}\, s(x;\theta) = -\frac{x}{\theta^2} - \frac{1-x}{(1-\theta)^2}. \]
Taking \(-\mathbb{E}[\,\cdot\,]\) with \(\mathbb{E}[x]=\theta\) collapses the two fractions to \(\tfrac{1}{\theta} + \tfrac{1}{1-\theta}\), so the per-flip Fisher information and its inverse are
\[ I(\theta) = \frac{1}{\theta} + \frac{1}{1-\theta} = \frac{1}{\theta(1-\theta)}, \qquad \frac{I(\theta)^{-1}}{n} = \frac{\theta(1-\theta)}{n}. \tag{26.3.8}\]
Note that the per-flip information \(1/(\theta(1-\theta))\) is smallest at \(\theta = 1/2\): fair coins are the hardest to pin down, each flip carrying the least information about the bias. The prediction of Equation 26.3.5 is that the MLE \(\hat\theta = n_H/n\) over \(n\) flips has variance about \(\theta(1-\theta)/n\), which, by Equation 26.3.8, is also the Cramér–Rao floor Equation 26.3.7 for this problem. We can check that this is exactly right here (\(\hat\theta\) is a scaled binomial count, whose variance is \(\theta(1-\theta)/n\) on the nose) and confirm it by simulation: run the \(n\)-flip experiment many times and compare the spread of the resulting estimates to the floor.
# NumPy simulation of the coin MLE's variance.
theta_star, n, trials = 0.3, 200, 20000
rng = onp.random.default_rng(0)
# Each trial: flip the coin n times; the MLE is the observed head fraction.
heads = rng.binomial(n, theta_star, size=trials)
theta_hats = heads / n
empirical_var = theta_hats.var()
cramer_rao = theta_star * (1 - theta_star) / n # = I(theta)^-1 / n
empirical_var, cramer_rao(np.float64(0.0010385371777499998), 0.00105)(np.float64(0.0010385371777499998), 0.00105)(np.float64(0.0010385371777499998), 0.00105)(0.00103853717775, 0.00105)The empirical variance of the maximum-likelihood estimates matches the Cramér–Rao floor \(\theta(1-\theta)/n\) to two digits: the coin’s MLE is as tight as Equation 26.3.7 permits. This is asymptotic efficiency made concrete, and it is why the curvature of a loss surface (its Fisher information, or empirically its Hessian) tells us how trustworthy a fitted parameter is.
26.3.4 MAP Estimation: Priors as Regularizers
We dropped the prior \(P(\boldsymbol{\theta})\) by declaring it flat. Keeping it instead turns maximum likelihood into maximum a posteriori (MAP) estimation, and the prior reappears as a regularizer. Maximizing the full posterior Equation 26.3.9 and taking negative logs,
\[ \hat{\boldsymbol{\theta}}_{\textrm{MAP}} = \mathop{\mathrm{argmax}}_{\boldsymbol{\theta}}\, P(X\mid\boldsymbol{\theta})\,P(\boldsymbol{\theta}) = \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}}\, \Big[\underbrace{-\textstyle\sum_i \log p(x_i\mid\boldsymbol{\theta})}_{\textrm{data fit (NLL)}} \;\underbrace{-\log P(\boldsymbol{\theta})}_{\textrm{regularizer}}\Big]. \tag{26.3.9}\]
The objective is the familiar data-fit term plus a penalty: the negative log-prior. A prior concentrated near small parameters penalizes large ones, which is precisely what a regularizer does.
26.3.4.1 Gaussian Priors Are Weight Decay
Proposition (Gaussian prior = \(L_2\) / weight decay). A Gaussian prior \(\boldsymbol{\theta}\sim\mathcal{N}(\mathbf 0,\tau^2 \mathbf I)\) contributes the penalty \(\tfrac{1}{2\tau^2}\lVert\boldsymbol{\theta}\rVert_2^2\) to the MAP objective. With Gaussian-noise data (variance \(\sigma^2\)), MAP estimation is \(L_2\)-regularized least squares with weight-decay strength \(\lambda = \sigma^2/\tau^2\).
Proof. The density of \(\mathcal{N}(\mathbf 0,\tau^2\mathbf I)\) is proportional to \(\exp\!\big(-\lVert\boldsymbol{\theta}\rVert_2^2/(2\tau^2)\big)\), so its negative log-prior is \(\tfrac{1}{2\tau^2}\lVert\boldsymbol{\theta}\rVert_2^2\) plus a constant. Insert this and the Gaussian NLL Equation 26.3.4 into the MAP objective Equation 26.3.9:
\[ \hat{\boldsymbol{\theta}}_{\textrm{MAP}} = \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}}\; \frac{1}{2\sigma^2}\sum_{i}(y_i-\hat y_i)^2 + \frac{1}{2\tau^2}\lVert\boldsymbol{\theta}\rVert_2^2 . \]
Multiplying by the positive constant \(2\sigma^2\) leaves \(\sum_i(y_i-\hat y_i)^2 + \tfrac{\sigma^2}{\tau^2}\lVert\boldsymbol{\theta}\rVert_2^2\), which is ridge regression (Section 2.7) with \(\lambda=\sigma^2/\tau^2\). \(\blacksquare\)
This is the probabilistic origin of weight decay: the regularizer we add by hand in Section 2.7, and read as a norm constraint in Section 25.4.4.5, is the negative log of a Gaussian belief that weights should be small. The same construction with a Laplace prior \(p(\boldsymbol{\theta})\propto e^{-\lVert\boldsymbol{\theta}\rVert_1/b}\) gives the \(L_1\) penalty \(\tfrac{1}{b}\lVert\boldsymbol{\theta}\rVert_1\), whose kink at the origin produces exact zeros, i.e. sparsity. Figure 26.3.3 shows the tug-of-war: the NLL pulls toward the MLE, the log-prior pulls toward the prior mean, and the MAP estimate sits between them.
The figure also reads off the two limits. As the prior flattens (\(\tau\to\infty\), or \(b\to\infty\)) the penalty vanishes and \(\hat{\boldsymbol{\theta}}_{\textrm{MAP}} \to \hat{\boldsymbol{\theta}}_{\textrm{MLE}}\), the flat-prior assumption we started with. And as the data grows the NLL term deepens faster than the fixed prior, so the data term dominates and MAP again approaches the MLE. A bookkeeping footnote makes the second limit quantitative: deep-learning code minimizes the mean per-example loss rather than the sum, and dividing the MAP objective by \(n\) rescales the weight decay to \(\lambda = \sigma^2/(n\tau^2)\): the prior’s pull on the averaged loss literally fades as \(1/n\). Regularization matters most precisely when data is scarce.
26.3.4.2 A Beta Prior on the Coin
Weight decay is MAP for regression; we now give the running coin its MAP treatment, completing the example. A prior on a coin’s bias must live on \([0,1]\), and the natural family is the Beta distribution \(p(\theta)\propto\theta^{a-1}(1-\theta)^{b-1}\) with shape parameters \(a, b > 0\), the conjugate prior of Section 26.2: multiplying by the likelihood Equation 26.3.2 just shifts the exponents by the observed counts, so the posterior is again a Beta (Equation 26.2.18),
\[ P(\theta\mid X) \;\propto\; \theta^{\,n_H + a - 1}\,(1-\theta)^{\,n_T + b - 1}. \]
MAP asks for the mode of this density. That is the same single-variable maximization we solved for the likelihood itself: the posterior has the shape \(\theta^{\alpha}(1-\theta)^{\beta}\) with \(\alpha = n_H+a-1\) and \(\beta = n_T+b-1\), and whenever both exponents are positive (e.g. \(a, b \ge 1\) on our data) the maximum is interior. The answer is therefore the same ratio, “head exponent over total exponent”:
\[ \hat\theta_{\textrm{MAP}} \;=\; \frac{n_H + a - 1}{n + a + b - 2}, \qquad n = n_H + n_T. \tag{26.3.10}\]
The reading is pseudo-counts: the prior behaves exactly like \(a-1\) phantom heads and \(b-1\) phantom tails appended to the data before taking the empirical frequency. A \(\mathrm{Beta}(2,2)\) prior adds one phantom flip of each kind, turning our \(9\)-of-\(13\) coin into \((9+1)/(13+2) = 2/3\), pulled from \(9/13\approx 0.692\) toward \(1/2\); a \(\mathrm{Beta}(30,30)\) prior would all but pin the estimate near a fair coin until the data outweighed its \(58\) phantom flips. And as \(n\) grows, any fixed stock of phantom flips washes out: \(\hat\theta_{\textrm{MAP}}\to\hat\theta_{\textrm{MLE}}\), the data-swamps-the-prior limit we just met.
Two subtleties in Equation 26.3.10 are easy to get wrong. First, the flat prior \(\mathrm{Beta}(1,1)\), uniform on \([0,1]\), gives \(\hat\theta_{\textrm{MAP}} = n_H/n\): MAP under a uniform prior is exactly the MLE, zero phantom flips, no smoothing at all. In particular the “add-one” (Laplace) smoothing rule \((n_H+1)/(n+2)\) is not the uniform-prior MAP; read as pseudo-counts, it is the MAP under \(\mathrm{Beta}(2,2)\).
Second, the mode is not the mean. The posterior mean is
\[ \mathbb{E}[\theta\mid X] \;=\; \frac{n_H + a}{n + a + b}, \]
which under the flat prior becomes \((n_H+1)/(n+2)\): Laplace’s rule of succession (Laplace 1814), precisely the add-one rule. So add-one smoothing is a Bayes estimate under the uniform prior, just a different summary: the posterior mean, not the posterior mode. (Numerically it coincides with the \(\mathrm{Beta}(2,2)\) mode, which is why the two are so often conflated.) Naive Bayes (Section 26.7) leans on exactly this smoothing to keep never-observed feature–class pairs from annihilating its products; the pseudo-count accounting built here is what justifies it.
26.3.4.3 The Posterior Mode Is Not the Posterior
One caution before we move on: MAP is not the same as “being Bayesian.” MAP restores the prior but still reports a single point: the mode of the posterior \(P(\boldsymbol{\theta}\mid X)\), the same \(\mathop{\mathrm{argmax}}\) as in Equation 26.3.9. Genuine Bayesian inference keeps the entire posterior distribution and integrates over it, predicting with the posterior-averaged \(\int p(x\mid\boldsymbol{\theta})\,P(\boldsymbol{\theta}\mid X)\,d\boldsymbol{\theta}\) rather than plugging in one \(\hat{\boldsymbol{\theta}}_{\textrm{MAP}}\). That average propagates parameter uncertainty into the prediction, where the mode discards it. The mode is not even reparameterization-invariant, since a nonlinear change of variables moves the peak of a density but not its integral. Maximum likelihood does not suffer this wobble: it is equivariant under reparameterization, meaning that if \(\hat\theta\) is the MLE of \(\theta\), then \(g(\hat\theta)\) is the MLE of \(g(\theta)\) for any one-to-one \(g\), so the MLE of the coin’s log-odds is simply \(\log\frac{9/13}{4/13} = \log\frac{9}{4}\) (Wasserman 2013). The reason is that relabeling the parameter axis moves no maximum: likelihood values are attached to points, while densities carry a Jacobian. MAP inherits the prior density’s Jacobian, and with it the dependence on how the parameter happens to be written down. Marginalizing over the posterior is its own subject, beyond our scope here (Murphy 2022; MacKay 2003); the point is only that MAP is one more point estimate (maximum likelihood with a penalty) and should not be mistaken for the full Bayesian treatment.
26.3.5 Latent Variables, EM, and the ELBO
Every likelihood so far has been a product of terms we could log and differentiate directly. The most interesting generative models are not so obliging. They posit a latent variable \(z\) (an unobserved cluster identity, a topic, the “seed” a generator decodes into an image) and define the joint distribution
\[ p(x, z; \boldsymbol{\theta}) = p(z; \boldsymbol{\theta})\, p(x \mid z; \boldsymbol{\theta}), \]
of which we observe only \(x\). The likelihood of an observation is then a marginal, obtained by summing (for continuous latents, integrating) the latent out:
\[ p(x;\boldsymbol{\theta}) = \sum_{z} p(z;\boldsymbol{\theta})\, p(x\mid z;\boldsymbol{\theta}). \tag{26.3.11}\]
26.3.5.1 Why Latent Variables Break the Recipe
The canonical example is the Gaussian mixture model (GMM): the latent \(z \in \{1,\ldots,K\}\) picks a component with probability \(P(z = k) = \pi_k\), and the observation is then drawn from that component, \(x \mid z{=}k \sim \mathcal{N}(\mu_k, \sigma_k^2)\). The data’s NLL is
\[ \ell(\boldsymbol{\theta}) = -\sum_{i=1}^n \log\Big(\sum_{k=1}^K \pi_k\, \mathcal{N}(x_i;\, \mu_k, \sigma_k^2)\Big), \qquad \boldsymbol{\theta} = (\pi_k, \mu_k, \sigma_k)_{k=1}^K. \]
The trouble is the sum inside the logarithm. Until now every log reached an individual density and the NLL split into per-example, per-parameter pieces; here the log stops at the mixture, every parameter appears inside every term through the shared sum, and no closed form solves the first-order conditions. Contrast the complete-data problem: if an oracle revealed each \(z_i\), the likelihood would factor by component, and maximum likelihood would collapse to the estimators we already own: the fraction of points in component \(k\) for \(\pi_k\) (the coin), and the per-component sample mean and variance for \(\mu_k, \sigma_k^2\) (the Gaussian). Plain gradient descent on \(\ell(\boldsymbol{\theta})\) remains possible, but the structure of the problem (easy if complete, hard only because \(z\) is missing) suggests something better: guess the missing labels probabilistically, then solve the easy problem.
26.3.5.2 The Evidence Lower Bound
Let \(q(z)\) be any distribution over the latent with \(q(z) > 0\) wherever \(p(x,z;\boldsymbol{\theta}) > 0\), our probabilistic guess. Insert it into the marginal and apply Jensen’s inequality (\(\log\) is concave, so \(\log \mathbb{E}[\cdot] \ge \mathbb{E}[\log(\cdot)]\); Section 25.3.3):
\[ \log p(x;\boldsymbol{\theta}) = \log \sum_z q(z)\, \frac{p(x,z;\boldsymbol{\theta})}{q(z)} \;\ge\; \sum_z q(z) \log \frac{p(x,z;\boldsymbol{\theta})}{q(z)} \;=:\; \mathcal{L}(q, \boldsymbol{\theta}). \tag{26.3.12}\]
The right-hand side is the evidence lower bound ELBO. The name comes from calling the marginal log-likelihood \(\log p(x;\boldsymbol{\theta})\) the evidence: it plays for the latent \(z\) exactly the role the denominator \(P(X)\) played for the parameters in Bayes’ rule at the start of this section: the marginal probability of what we actually observed. Expanding the fraction splits the bound into two readable pieces,
\[ \mathcal{L}(q,\boldsymbol{\theta}) = \mathbb{E}_{q}\big[\log p(x,z;\boldsymbol{\theta})\big] + H(q), \]
the expected complete-data log-likelihood (the easy objective, averaged over the guess) plus the entropy of the guess, which does not involve \(\boldsymbol{\theta}\) at all (Section 27.1). How much does the bound give away? Exactly a KL divergence:
Proposition (the ELBO gap is a KL). For every \(q\) and every \(\boldsymbol{\theta}\),
\[ \log p(x;\boldsymbol{\theta}) = \mathcal{L}(q,\boldsymbol{\theta}) + D_{\textrm{KL}}\big(q(z)\,\|\,p(z\mid x;\boldsymbol{\theta})\big). \tag{26.3.13}\]
Hence \(\mathcal{L}(q,\boldsymbol{\theta}) \le \log p(x;\boldsymbol{\theta})\), with equality if and only if \(q\) is the posterior \(p(z\mid x;\boldsymbol{\theta})\).
Proof. Expand the KL term using Bayes’ rule, \(p(z\mid x;\boldsymbol{\theta}) = p(x,z;\boldsymbol{\theta})/p(x;\boldsymbol{\theta})\):
\[ D_{\textrm{KL}}\big(q \,\|\, p(z\mid x;\boldsymbol{\theta})\big) = \mathbb{E}_q\!\left[\log \frac{q(z)\, p(x;\boldsymbol{\theta})}{p(x,z;\boldsymbol{\theta})}\right] = \log p(x;\boldsymbol{\theta}) - \mathcal{L}(q,\boldsymbol{\theta}), \]
where \(\log p(x;\boldsymbol{\theta})\) exits the expectation because it does not involve \(z\). Rearranging gives Equation 26.3.13; since the KL divergence is nonnegative and vanishes exactly when its arguments agree, the bound and its equality condition follow. \(\blacksquare\)
Figure 26.3.4 draws the resulting geometry: for a fixed guess the ELBO is a curve running everywhere below the evidence, and choosing \(q\) well pinches the two together at the current parameters.
26.3.5.3 Expectation–Maximization
The proposition turns the intractable problem into coordinate ascent on \(\mathcal{L}(q, \boldsymbol{\theta})\): maximize over one of the two arguments while the other is held fixed, then swap, and repeat. For a dataset \(X = \{x_1,\ldots,x_n\}\) of independent observations the bound applies example by example and sums, with one guess \(q_i\) per data point. The expectation–maximization (EM) algorithm alternates the two coordinates (Dempster et al. 1977):
- E-step. Hold \(\boldsymbol{\theta}^{(t)}\) fixed and maximize over the guesses: by the equality condition of Equation 26.3.13, the best \(q_i\) is the posterior, \(q_i^{(t)}(z) = p(z \mid x_i; \boldsymbol{\theta}^{(t)})\). This closes every gap, making the bound tight: \(\sum_i \mathcal{L}(q_i^{(t)}, \boldsymbol{\theta}^{(t)}) = \log p(X;\boldsymbol{\theta}^{(t)})\).
- M-step. Hold the guesses fixed and maximize over the parameters: \(\boldsymbol{\theta}^{(t+1)} = \mathop{\mathrm{argmax}}_{\boldsymbol{\theta}} \sum_i \mathbb{E}_{q_i^{(t)}}\big[\log p(x_i, z; \boldsymbol{\theta})\big]\) (the entropy terms do not involve \(\boldsymbol{\theta}\)). This is the easy, complete-data objective, with soft assignments standing in for the oracle’s labels.
Proposition (EM never decreases the likelihood). Each EM iteration satisfies \(\log p(X;\boldsymbol{\theta}^{(t+1)}) \ge \log p(X;\boldsymbol{\theta}^{(t)})\).
Proof. Chain three facts:
\[ \log p(X;\boldsymbol{\theta}^{(t+1)}) \;\ge\; \sum_i \mathcal{L}\big(q_i^{(t)}, \boldsymbol{\theta}^{(t+1)}\big) \;\ge\; \sum_i \mathcal{L}\big(q_i^{(t)}, \boldsymbol{\theta}^{(t)}\big) \;=\; \log p(X;\boldsymbol{\theta}^{(t)}). \]
The first inequality is the ELBO bound Equation 26.3.12 at the new parameters; the second is the definition of the M-step; the equality is the tightness the E-step bought. \(\blacksquare\)
For the Gaussian mixture both steps are closed-form. The E-step computes responsibilities, the posterior probability that point \(i\) came from component \(k\),
\[ r_{ik} = \frac{\pi_k\, \mathcal{N}(x_i;\, \mu_k, \sigma_k^2)} {\sum_{j} \pi_j\, \mathcal{N}(x_i;\, \mu_j, \sigma_j^2)}, \]
and the M-step refits each component by responsibility-weighted counting and averaging: with effective counts \(n_k = \sum_i r_{ik}\),
\[ \pi_k = \frac{n_k}{n}, \qquad \mu_k = \frac{1}{n_k}\sum_{i} r_{ik}\, x_i, \qquad \sigma_k^2 = \frac{1}{n_k}\sum_{i} r_{ik}\,(x_i - \mu_k)^2. \tag{26.3.14}\]
These are the coin’s counting estimator and the Gaussian’s mean and variance estimators, softened: instead of voting for one component, each point spreads its single vote across components in proportion to its responsibilities (exercise 10 derives the \(\mu_k\) update). The cell below runs EM on a two-component mixture in plain NumPy and prints the log-likelihood as it climbs.
# EM for a two-component 1-D Gaussian mixture, in plain numpy.
rng = onp.random.default_rng(42)
# 300 points from N(-2, 0.8^2) and 200 from N(1.5, 0.6^2)
x = onp.concatenate([rng.normal(-2.0, 0.8, 300), rng.normal(1.5, 0.6, 200)])
def gauss(x, mu, sigma):
return onp.exp(-0.5 * ((x - mu) / sigma)**2) / (sigma * onp.sqrt(2 * onp.pi))
# Deliberately poor start: equal weights, overlapping components
pi, mu, sigma = onp.array([0.5, 0.5]), onp.array([-0.5, 0.5]), onp.array([1.0, 1.0])
for step in range(31):
# E-step: responsibilities r[i, k], the posterior of component k given x_i
joint = pi * gauss(x[:, None], mu, sigma)
r = joint / joint.sum(axis=1, keepdims=True)
if step % 5 == 0: # the log-likelihood never decreases
print(f'step {step:2d}, log-likelihood {onp.log(joint.sum(axis=1)).sum():.2f}')
# M-step: refit each component by responsibility-weighted averaging
n_k = r.sum(axis=0)
pi, mu = n_k / len(x), (r * x[:, None]).sum(axis=0) / n_k
sigma = onp.sqrt((r * (x[:, None] - mu)**2).sum(axis=0) / n_k)
print(f'pi = {pi.round(3)}, mu = {mu.round(3)}, sigma = {sigma.round(3)}')step 0, log-likelihood -1289.53
step 5, log-likelihood -845.69
step 10, log-likelihood -845.45
step 15, log-likelihood -845.45
step 20, log-likelihood -845.45
step 25, log-likelihood -845.45
step 30, log-likelihood -845.45
pi = [0.594 0.406], mu = [-2.055 1.494], sigma = [0.715 0.625]The trace behaves exactly as the proposition demands: the log-likelihood rises monotonically, gains most of its ground in the first few iterations, and then flatlines at a stationary point; the recovered weights, means, and scales land near the generating values \((0.6, 0.4)\), \((-2, 1.5)\), \((0.8, 0.6)\), the small residual gaps being finite-sample error, not a failure of EM. Two cautions temper the success. EM is ascent, not global optimization: the mixture likelihood is multimodal (at minimum because relabeling the components permutes the parameters without changing the model, the same identifiability failure we met in the estimator theory above), and EM can stall at a poor local maximum, so practice restarts it from several initializations. And strictly speaking the GMM likelihood is unbounded: let one component collapse onto a single data point with \(\sigma_k \to 0\). Monotonicity itself guarantees only nondecrease; because this objective is unbounded, the values may diverge rather than approach a finite limit. With a variance floor or prior making the relevant level set bounded, plus standard continuity and regularity conditions, EM limit points are stationary (Wu 1983).
The bound, not the closed-form M-step, is the part that scales. A variational autoencoder (Kingma and Welling 2014) keeps the ELBO but replaces the per-point E-step with a neural network \(q_{\boldsymbol{\phi}}(z\mid x)\), an amortized guess trained jointly with the generator by stochastic gradient ascent on the ELBO, because for a deep decoder the true posterior is no longer available in closed form. Diffusion models walk further down the same road: their training objective is the ELBO of a deep chain of latent variables (Section 28.4). And the two divergences in this section are deliberately mirrored: Section 27.2 recasts maximum likelihood itself as forward-KL minimization, while the ELBO gap Equation 26.3.13 measures the guess against the posterior in the reverse orientation, \(D_{\textrm{KL}}(q\,\|\,p)\), a distinction with real modelling consequences, taken up in detail there.
26.3.6 Summary
- The maximum likelihood principle directly picks the parameter that maximizes the observed-data probability or density, \(\mathop{\mathrm{argmax}}_{\boldsymbol{\theta}}P(X\mid\boldsymbol{\theta})\). It does not require a prior; MAP is a separate posterior-mode rule that coincides with MLE under a constant prior in the chosen parameterization.
- We optimize the negative log-likelihood \(-\sum_i \log p(x_i\mid\boldsymbol{\theta})\): the log fixes underflow and makes gradients an additive sum, while the sign flip turns “maximize probability” into “minimize a loss.”
- For discrete outcomes, the average NLL is the cross-entropy from the empirical frequencies to the model and differs from the corresponding KL by a fixed empirical entropy. For continuous data, use the average negative log-density at the sample; the clean KL-projection statement is a population statement under suitable integrability conditions. Categorical NLL is softmax cross-entropy, fixed-variance Gaussian NLL is mean squared error, and a Laplace model gives mean absolute error: picking a loss is picking a noise model.
- The MLE is consistent (given identifiability and uniform convergence) and, under the regularity conditions granted in the text, asymptotically normal with limiting variance the inverse Fisher information \(I(\boldsymbol{\theta}^\star)^{-1}/n\), so it attains the Cramér–Rao bound Equation 26.3.7 and is asymptotically efficient. At finite \(n\) it is generally only asymptotically unbiased (Section 26.5).
- MAP estimation restores the prior, adding its negative log as a penalty: a Gaussian prior is exactly \(L_2\) / weight decay (\(\lambda=\sigma^2/\tau^2\); for the averaged loss, \(\sigma^2/(n\tau^2)\)), a Laplace prior is \(L_1\) / sparsity, and a Beta prior on the coin adds pseudo-counts, with MAP under the flat Beta(1,1) prior being the plain MLE, while add-one smoothing is the posterior mean (rule of succession). MAP \(\to\) MLE as the prior flattens or the data grows, and it is still a point estimate (the posterior mode), not full Bayesian inference, which integrates over the posterior.
- When the model has latent variables the likelihood is a marginal with a sum inside the log. The ELBO lower-bounds it for every guess \(q\), with gap exactly \(D_{\textrm{KL}}(q \,\|\, p(z\mid x;\boldsymbol{\theta}))\); EM alternates closing the gap (E-step: \(q\) = posterior) with climbing the bound (M-step), and never decreases the likelihood. VAEs and diffusion models train on the same bound with amortized guesses.
- Maximum likelihood extends to continuous variables unchanged: probabilities become densities, and the matching-tolerance \(\epsilon\) drops out as a constant.
26.3.7 Exercises
- A non-negative random variable has density \(\alpha e^{-\alpha x}\) for some \(\alpha>0\). From the single observation \(x=3\), find the maximum likelihood estimate of \(\alpha\). Generalize to a sample \(\{x_i\}_{i=1}^N\) and show \(\hat\alpha = 1/\bar x\).
- Given a sample \(\{x_i\}_{i=1}^N\) from a Gaussian with unknown mean and variance \(1\), show the maximum likelihood estimate of the mean is the sample average \(\bar x\). (Hint: minimize the Gaussian NLL of Section 26.3.2.3.)
- Show directly that a fixed-variance Gaussian NLL equals the mean squared error up to an additive constant, and identify the constant.
- For a \(K\)-class classifier with one-hot label \(\mathbf y\) and softmax prediction \(\hat{\mathbf p} = \mathrm{softmax}(\mathbf o)\), where the vector \(\mathbf o\) collects the logits, the raw scores fed to the softmax (Section 3.1), show the categorical NLL \(-\sum_k y_k\log\hat p_k\) is the cross-entropy, and that its gradient with respect to the logits is \(\hat{\mathbf p} - \mathbf y\) (compare Equation 24.3.7).
- Show that a Gaussian prior on the regression weights recovers ridge regression, and read off the weight-decay strength in terms of the noise and prior variances.
- Show that MAP estimation with a Laplace prior gives an \(L_1\) penalty, and argue why its kink at the origin can force coefficients to be exactly zero.
- Prove that \(\hat{\boldsymbol{\theta}}_{\textrm{MAP}} \to \hat{\boldsymbol{\theta}}_{\textrm{MLE}}\) as the prior variance \(\tau^2\to\infty\), and explain why the same happens as the dataset size \(n\to\infty\) for fixed \(\tau\).
- For a single Gaussian observation with known variance \(\sigma^2\) and unknown mean \(\mu\), compute the score \(\frac{d}{d\mu}\log p(x\mid\mu)\) and show the Fisher information is \(I(\mu)=1/\sigma^2\). Conclude that the sample-mean MLE over \(n\) observations has variance \(\sigma^2/n\), attaining the Cramér–Rao bound Equation 26.3.7 exactly. (Hint: use Equation 26.3.6.)
- Extend the asymptotic-normality demonstration in the text. First rerun it with \(n = 10\) in place of \(400\): how many distinct values can \(\hat\theta\) take, and what does the “histogram” become? Then restore \(n = 400\) but move \(\theta^\star\) to \(0.95\) and compare the match to the predicted Gaussian. Explain why convergence is slower near the boundary, considering both the skew of the binomial there and the shrinking variance \(\theta^\star(1-\theta^\star)\).
- Derive the GMM M-step for the means: holding the responsibilities \(r_{ik}\) fixed, show that maximizing the expected complete-data log-likelihood \(\sum_i \sum_k r_{ik}\big[\log \pi_k + \log \mathcal{N}(x_i;\,\mu_k,\sigma_k^2)\big]\) over \(\mu_k\) yields the responsibility-weighted mean of
- Then derive the \(\pi_k\) update by maximizing over the weights subject to \(\sum_k \pi_k = 1\). (Hint: use a Lagrange multiplier for the constraint.)