Maximum likelihood, the exponential family and properties of estimators

Before this lesson: lesson 03 (Bernoulli, Poisson, exponential, Gaussian and gamma distributions), lesson 04 (covariance matrices are positive semidefinite; the sample variance $S^2$) and lesson 06 (the law of large numbers). From other tracks: gradients and Hessians (calculus 02, calculus 03); the KL divergence (information theory 03) is defined again in section 3. Time: about 75 minutes, plus the notebook.

---

Why this matters #

• Training losses are likelihoods. The mean squared error is the negative log-likelihood (NLL) of a Gaussian model. Binary and softmax cross-entropy are the NLLs of Bernoulli and categorical models. Minimising them is maximum likelihood (information theory 04).
• Error bars and comparisons. Bias, variance and MSE are the language for comparing estimators: why the sample variance divides by $n - 1$, why a regularised (shrunk) estimate can beat the MLE, why a model evaluated with few random seeds looks more stable than it is.
• Generalised linear models. Linear, logistic and Poisson regression are maximum likelihood in an exponential family. This is why their losses are convex and why their gradients all look like "prediction minus target".
• Science. Fitting a gamma distribution to rainfall amounts, an exponential distribution to the waiting times between events, or a Poisson rate to counts of photons or earthquakes is maximum likelihood.

1. Predict first #

2. The picture #

The likelihood of a coin #

You flip a coin $n = 10$ times and see $k = 7$ heads. Model: each flip is $x_i \sim \text{Bernoulli}(p)$, independent, with $x_i = 1$ for heads. The probability of the observed sequence is

Read this formula in two ways.

• As a probability: fix $p$ and vary the data. For $p = 0.7$, it gives the probability of each possible sequence of 10 flips. These numbers sum to 1 over all $2^{10}$ sequences.
• As a likelihood: fix the data (7 heads) and vary $p$. This gives a curve over $p \in [0, 1]$:

The likelihood $L(\theta)$ is the probability (or density) of the observed data, seen as a function of the parameter $\theta$. It is the same formula, with a different variable. The maximum likelihood estimate is the value of $\theta$ where $L$ is largest.

Some values: $L(0.7) = 0.7^7 \cdot 0.3^3 \approx 0.00222$ and $L(0.5) = 0.5^{10} \approx 0.00098$. So the data are about 2.3 times more probable under $p = 0.7$ than under a fair coin. The number $0.00222$ itself has no meaning on its own; only ratios of likelihoods do.

Interactive figure: loading…

Why we take logs. The likelihood of i.i.d. data is a product of $n$ numbers, often much smaller than 1. For 1000 flips with 700 heads, $L(0.7)$ is about $10^{-265}$; with 2000 flips the product underflows to 0 in double precision. The log-likelihood $\ell(\theta) = \log L(\theta)$ turns the product into a sum,

and has the same maximiser, because $\log$ is increasing. Set the derivative to zero:

Check that it is a maximum: $\ell''(p) = -\frac{k}{p^2} - \frac{n-k}{(1-p)^2} < 0$ for all $p \in (0, 1)$. At $\hat p = 0.7$, $\ell''(\hat p) = -\frac{7}{0.49} - \frac{3}{0.09} \approx -47.6$, which equals $-\frac{n}{\hat p(1 - \hat p)} = -\frac{10}{0.21}$. The size of this second derivative is the curvature of the peak. It grows in proportion to $n$.

The likelihood is not a density in $p$. The area under $L(p)$ is $\int_0^1 p^7(1-p)^3\,dp = \frac{7!\,3!}{11!} = \frac{1}{1320}$, not 1. To get a probability distribution over $p$ you need a prior and Bayes' rule; that is lesson 10.

Fitting a Gaussian: every point votes #

Now the data are real numbers $x_1, \dots, x_n$ and the model is $\mathcal N(\mu, \sigma^2)$. Each point contributes its log-density:

Think of $-\log p(x_i)$ as the "cost" of point $i$. A point far from $\mu$ (measured in units of $\sigma$) costs a lot. A large $\sigma$ makes far points cheap, but it makes every point pay the term $\log\sigma$. The MLE balances the two.

Interactive figure: loading…

The minimum of the total cost is at $\hat\mu = \bar x$ and $\hat\sigma^2 = \frac{1}{n}\sum_i (x_i - \bar x)^2$ (section 3 proves it). In words: the fitted Gaussian has the same mean and the same average squared deviation as the data.

3. Definitions and main results #

Likelihood and the maximum likelihood estimate #

Setting. A model is a family of PMFs or densities $p(x \mid \theta)$, indexed by a parameter $\theta \in \Theta$ (a number or a vector). The data $x_1, \dots, x_n$ are i.i.d. from $p(x \mid \theta)$ for one true value of $\theta$.

• The likelihood is $L(\theta) = \prod_{i=1}^n p(x_i \mid \theta)$, a function of $\theta$ for fixed data.
• The log-likelihood is $\ell(\theta) = \sum_{i=1}^n \log p(x_i \mid \theta)$.
• The maximum likelihood estimate is $\hat\theta_{\text{MLE}} = \arg\max_{\theta \in \Theta} \ell(\theta)$, when the maximum exists.

In practice, set the gradient $\nabla\ell(\theta)$ to zero (the likelihood equations) and check that the stationary point is a maximum: a negative second derivative, or a negative definite Hessian (calculus 03). Write $\bar x = \frac{1}{n}\sum_i x_i$.

Data	$\ell(\theta)$, up to constants	Solve $\ell'(\theta) = 0$	MLE
Bernoulli($p$)	$n\bar x\log p + n(1-\bar x)\log(1-p)$	$\frac{n\bar x}{p} = \frac{n(1-\bar x)}{1-p}$	$\hat p = \bar x$
Poisson($\lambda$)	$n\bar x\log\lambda - n\lambda$	$\frac{n\bar x}{\lambda} - n = 0$	$\hat\lambda = \bar x$
Exponential($\lambda$), rate $\lambda$	$n\log\lambda - \lambda n\bar x$	$\frac{n}{\lambda} - n\bar x = 0$	$\hat\lambda = 1/\bar x$
$\mathcal N(\mu, v)$ with $v = \sigma^2$	$-\frac{n}{2}\log v - \frac{1}{2v}\sum_i (x_i - \mu)^2$	see below	$\hat\mu = \bar x$, $\hat v = \frac{1}{n}\sum_i (x_i - \bar x)^2$

The second derivatives are negative in all four cases: $-\frac{n\bar x}{\lambda^2}$ (Poisson, when $\bar x > 0$), $-\frac{n}{\lambda^2}$ (exponential), and $-\frac{n\bar x}{p^2} - \frac{n(1-\bar x)}{(1-p)^2}$ (Bernoulli).

Gaussian, step by step. We maximise over $\mu$ and over the variance $v = \sigma^2$ (not $\sigma$; the answer is the same by invariance, below).

1. $\partial_\mu \ell = \frac{1}{v}\sum_i (x_i - \mu)$. It is zero when $\sum_i x_i = n\mu$, so $\hat\mu = \bar x$, for every $v$.
2. $\partial_v \ell = -\frac{n}{2v} + \frac{1}{2v^2}\sum_i (x_i - \hat\mu)^2$. Multiply by $2v^2$: $-nv + \sum_i (x_i - \bar x)^2 = 0$, so $\hat v = \frac{1}{n}\sum_i (x_i - \bar x)^2$.
3. It is a maximum: for fixed $v$, $\ell$ is a downward parabola in $\mu$. With $\mu = \bar x$ and $D = \sum_i (x_i - \bar x)^2 > 0$, the function $-\frac{n}{2}\log v - \frac{D}{2v}$ increases for $v < D/n$ and decreases after, so $v = D/n$ is its maximum. (If all $x_i$ are equal, $D = 0$ and $\ell \to \infty$ as $v \to 0$: no MLE exists.)

Invariance of the MLE. Let $\phi = g(\theta)$ for a one-to-one map $g$. The likelihood in the new parameter is $\tilde L(\phi) = L(g^{-1}(\phi))$. It takes the same values, only relabelled, so its maximum is at $\hat\phi = g(\hat\theta)$. Examples: the MLE of $\sigma$ is $\sqrt{\hat v}$; the MLE of the log-odds $\log\frac{p}{1-p}$ is $\log\frac{\hat p}{1-\hat p}$; the MLE of the mean waiting time $1/\lambda$ is $1/\hat\lambda = \bar x$.

Estimators are random variables #

An estimator is a rule $\hat\theta = \hat\theta(X_1, \dots, X_n)$ that turns data into a guess. Before the data are drawn, $\hat\theta$ is a random variable with its own distribution, the sampling distribution. For a true value $\theta$:

• the bias is $\operatorname{bias}(\hat\theta) = \mathbb{E}_\theta[\hat\theta] - \theta$; the estimator is unbiased if the bias is 0 for every $\theta$;
• the variance is $\operatorname{Var}_\theta(\hat\theta)$, and its square root is the standard error;
• the mean squared error is $\text{MSE}(\hat\theta) = \mathbb{E}_\theta[(\hat\theta - \theta)^2]$.

A dartboard is the usual picture: the bias is how far the centre of the darts is from the bull's-eye, the variance is how spread out the darts are, and the MSE is the average squared distance to the bull's-eye.

Theorem (bias–variance decomposition). $\text{MSE}(\hat\theta) = \operatorname{bias}(\hat\theta)^2 + \operatorname{Var}(\hat\theta)$.

Proof.

1. Write $m = \mathbb{E}[\hat\theta]$ and $b = m - \theta$ (a constant). Split the error: $\hat\theta - \theta = (\hat\theta - m) + b$.
2. Square: $(\hat\theta - \theta)^2 = (\hat\theta - m)^2 + 2b\,(\hat\theta - m) + b^2$.
3. Take expectations. $\mathbb{E}[\hat\theta - m] = 0$, so the middle term vanishes, and $\mathbb{E}[(\hat\theta - \theta)^2] = \operatorname{Var}(\hat\theta) + b^2$. $\blacksquare$

Example: three ways to estimate a variance. Let $X_1, \dots, X_n$ be i.i.d. $\mathcal N(\mu, \sigma^2)$ and $D = \sum_i (X_i - \bar X)^2$. Lesson 04 showed $\mathbb{E}[D] = (n-1)\sigma^2$. For Gaussian data, $D/\sigma^2$ has a chi-squared distribution with $n - 1$ degrees of freedom, so $\operatorname{Var}(D) = 2(n-1)\sigma^4$. Consider the estimators $D/c$:

Divide by	Name	Bias	Variance	MSE at $n = 10$, $\sigma^2 = 1$
$c = n - 1$	$S^2$, unbiased	$0$	$\frac{2\sigma^4}{n-1}$	$0.222$
$c = n$	$\hat\sigma^2_{\text{MLE}}$	$-\frac{\sigma^2}{n}$	$\frac{2(n-1)\sigma^4}{n^2}$	$0.190$
$c = n + 1$		$-\frac{2\sigma^2}{n+1}$	$\frac{2(n-1)\sigma^4}{(n+1)^2}$	$0.182$

The bias of the MLE comes from $\mathbb{E}[D/n] = \frac{n-1}{n}\sigma^2$. But the biased estimators have a smaller MSE: they pay a little bias for a larger cut in variance. For Gaussian data, $c = n + 1$ gives the smallest MSE of all estimators of the form $D/c$. Unbiased is not the same as best.

Interactive figure: loading…

Another bias: a nonlinear function. For exponential data, $\bar X$ is unbiased for the mean $1/\lambda$. But $\hat\lambda = 1/\bar X$ is not unbiased for $\lambda$. The sum $\sum_i X_i$ has a Gamma($n, \lambda$) distribution (lesson 03), and a short integral gives, for $n \ge 3$,

At $n = 10$ and $\lambda = 1$: mean $1.111$, variance $0.154$, MSE $0.167$. The bias is upward, as Jensen's inequality predicts: $t \mapsto 1/t$ is convex, so $\mathbb{E}[1/\bar X] \ge 1/\mathbb{E}[\bar X] = \lambda$.

Consistency. A sequence of estimators $\hat\theta_n$ (one for each sample size) is consistent if $\hat\theta_n \to \theta$ in probability for every $\theta$: $P(|\hat\theta_n - \theta| > \varepsilon) \to 0$ for every $\varepsilon > 0$.

• MSE → 0 is enough. By Markov's inequality applied to $(\hat\theta_n - \theta)^2$ (lesson 06), $P(|\hat\theta_n - \theta| > \varepsilon) \le \text{MSE}(\hat\theta_n)/\varepsilon^2 \to 0$.
• $\bar X$ is consistent (law of large numbers). $\hat\sigma^2_{\text{MLE}}$ is biased but consistent: its bias $-\sigma^2/n$ and its variance both go to 0. And $1/\bar X$ is consistent for $\lambda$, because $t \mapsto 1/t$ is continuous at $1/\lambda$ (continuous mapping).
• Unbiased does not imply consistent: $X_1$ alone is an unbiased estimate of the mean, but its spread never shrinks.
• Under regularity conditions, the MLE is consistent. The KL view below explains why.

Why maximum likelihood: the KL view #

The Kullback–Leibler divergence between distributions $p$ and $q$ is

with equality only if $p = q$ (Gibbs' inequality, a consequence of Jensen's inequality; information theory 03). Let $p_{\text{data}}$ be the true distribution of the data, and $p_\theta$ the model.

1. By the law of large numbers, $\frac{1}{n}\ell(\theta) = \frac{1}{n}\sum_i \log p_\theta(X_i) \to \mathbb{E}_{p_{\text{data}}}[\log p_\theta(X)]$.
2. Add and subtract $\log p_{\text{data}}$: $\mathbb{E}_{p_{\text{data}}}[\log p_\theta(X)] = \mathbb{E}_{p_{\text{data}}}[\log p_{\text{data}}(X)] - \mathrm{KL}(p_{\text{data}} \,\|\, p_\theta)$.
3. The first term does not depend on $\theta$. So maximising the expected log-likelihood is the same as minimising $\mathrm{KL}(p_{\text{data}} \,\|\, p_\theta)$, and the same as minimising the cross-entropy $-\mathbb{E}_{p_{\text{data}}}[\log p_\theta(X)]$.
4. For finite $n$ the same identity holds exactly with the empirical distribution $\hat p_n$ (mass $1/n$ on each data point) in place of $p_{\text{data}}$: the MLE minimises the cross-entropy between $\hat p_n$ and $p_\theta$.

Three consequences:

• Consistency. If the model is correct ($p_{\text{data}} = p_{\theta_0}$) and identifiable (different $\theta$ give different distributions), the KL is zero only at $\theta_0$. So the limit of $\frac{1}{n}\ell$ has its unique maximum at the truth. (Turning this into $\hat\theta_n \to \theta_0$ needs a uniform law of large numbers or concavity.)
• Misspecification. If no $\theta$ fits exactly, the MLE still converges, to the $\theta^*$ that minimises $\mathrm{KL}(p_{\text{data}} \,\|\, p_\theta)$: the closest model in this sense.
• Losses are NLLs. The average negative log-likelihood is the training loss:

Model for the target $y$ given input $x$	Negative log-likelihood of one example	Name of the loss
$y \sim \mathcal N(f(x), \sigma^2)$, $\sigma$ fixed	$\frac{(y - f(x))^2}{2\sigma^2} + \text{const}$	mean squared error
$y \sim \text{Bernoulli}(\sigma(z(x)))$, logit $z$	$-y\log\sigma(z) - (1-y)\log(1 - \sigma(z))$	binary cross-entropy
$y \sim \text{Categorical}(\operatorname{softmax}(z(x)))$	$-\log \operatorname{softmax}(z)_y$	softmax cross-entropy
$y \sim \mathcal N(\mu(x), \sigma(x)^2)$, both predicted	$\frac12\log\sigma^2 + \frac{(y - \mu)^2}{2\sigma^2} + \text{const}$	Gaussian NLL

The exponential family #

Start with the coin. Rewrite the Bernoulli PMF with $\exp$ and $\log$:

The log-probability is linear in the data $x$, with the coefficient $\eta = \log\frac{p}{1-p}$, the log-odds. The rest, $\log(1-p)$, does not depend on $x$. Solve for $p$: $p = \frac{e^\eta}{1 + e^\eta} = \sigma(\eta)$, the sigmoid (foundations 05). So $\log(1 - p) = -\log(1 + e^\eta)$, and

Differentiate $A$: $A'(\eta) = \frac{e^\eta}{1 + e^\eta} = p = \mathbb{E}[X]$, and $A''(\eta) = p(1-p) = \operatorname{Var}(X)$. The derivatives of the normaliser are the mean and the variance. This is not a coincidence.

Definition. A family of densities (or PMFs) is an exponential family if it can be written as

(a sum instead of the integral for discrete $x$). The parts:

• $\eta \in \mathbb{R}^d$ is the natural parameter;
• $T(x) \in \mathbb{R}^d$ is the sufficient statistic;
• $h(x) \ge 0$ is the base measure; it does not depend on $\eta$;
• $A(\eta)$ is the log-partition function; it makes the total mass 1.

Because $h$ does not depend on $\eta$, the support (where $p > 0$) is the same for every $\eta$.

Distribution	$\eta$	$T(x)$	$h(x)$	$A(\eta)$
Bernoulli($p$)	$\log\frac{p}{1-p}$	$x$	$1$	$\log(1 + e^\eta)$
Poisson($\lambda$)	$\log\lambda$	$x$	$1/x!$	$e^\eta$
Exponential($\lambda$)	$-\lambda$	$x$	$1$ on $x \ge 0$	$-\log(-\eta)$
$\mathcal N(\mu, \sigma^2)$	$\big(\frac{\mu}{\sigma^2},\, -\frac{1}{2\sigma^2}\big)$	$(x,\, x^2)$	$\frac{1}{\sqrt{2\pi}}$	$-\frac{\eta_1^2}{4\eta_2} - \frac12\log(-2\eta_2)$
Gamma($k, \beta$), rate $\beta$	$(k - 1,\, -\beta)$	$(\log x,\, x)$	$1$ on $x > 0$	$\log\Gamma(\eta_1 + 1) - (\eta_1 + 1)\log(-\eta_2)$

• Gaussian. Expand the square: $-\frac{(x-\mu)^2}{2\sigma^2} = \frac{\mu}{\sigma^2}x - \frac{1}{2\sigma^2}x^2 - \frac{\mu^2}{2\sigma^2}$. The normaliser is $\frac{\mu^2}{2\sigma^2} + \log\sigma$; written in $\eta$ it is the entry of the table.
• Gamma (rate form, lesson 03): $\frac{\beta^k x^{k-1} e^{-\beta x}}{\Gamma(k)} = \exp\big((k-1)\log x - \beta x - (\log\Gamma(k) - k\log\beta)\big)$. The exponential is the case $k = 1$.
• Not exponential families: Uniform($0, \theta$) (its support depends on $\theta$), the Student-t and the Cauchy with an unknown centre, and mixtures.

Theorem (the log-partition function gives the moments). For $\eta$ in the interior of the set where $A(\eta) < \infty$,

Proof: differentiate the normaliser

1. Write $Z(\eta) = e^{A(\eta)} = \int h(x)\, e^{\eta^\top T(x)}\, dx$. In the interior of the parameter set you may differentiate under the integral: $\nabla Z(\eta) = \int T(x)\, h(x)\, e^{\eta^\top T(x)}\, dx$.
2. $\nabla A = \nabla \log Z = \frac{\nabla Z}{Z} = \int T(x)\, h(x)\, e^{\eta^\top T(x) - A(\eta)}\, dx = \mathbb{E}_\eta[T(X)]$.
3. For the second derivative, note that $\nabla_\eta\, p(x \mid \eta) = p(x \mid \eta)\,\big(T(x) - \nabla A(\eta)\big)$. Differentiate step 2 once more:

$$\nabla^2 A = \int T(x)\,\big(T(x) - \nabla A(\eta)\big)^\top p(x \mid \eta)\, dx = \mathbb{E}[TT^\top] - \mathbb{E}[T]\,\mathbb{E}[T]^\top = \operatorname{Cov}_\eta[T(X)]. \quad\blacksquare$$

Consequences.

• $A$ is convex. A covariance matrix is positive semidefinite (lesson 04), so $\nabla^2 A \succeq 0$.
• The log-likelihood is concave in $\eta$. For i.i.d. data, $\ell(\eta) = \sum_i \log h(x_i) + \eta^\top \sum_i T(x_i) - nA(\eta)$: a linear function minus $n$ times a convex one. So any stationary point is a global maximum.
• The MLE matches moments. $\nabla \ell(\eta) = \sum_i T(x_i) - n\nabla A(\eta) = 0$ gives

The fitted model reproduces the sample averages of the sufficient statistic. For the Gaussian, $T = (x, x^2)$: the fitted model has the data's mean and mean square, which is exactly $\hat\mu = \bar x$, $\hat\sigma^2 = \overline{x^2} - \bar x^2$. All four rows of the MLE table above are this rule.

Gamma: matching needs a computer. For the gamma, $\nabla A = (\psi(k) - \log\beta,\ k/\beta) = (\mathbb{E}[\log X],\ \mathbb{E}[X])$, where $\psi = \Gamma'/\Gamma$ is the digamma function. So the MLE solves

Put $\hat\beta = \hat k/\bar x$ into the first equation: $\log\hat k - \psi(\hat k) = \log\bar x - \frac{1}{n}\sum_i\log x_i$. The right side is $\ge 0$ by Jensen's inequality ($\log$ is concave). There is no closed form; Newton's method solves it in a few steps (notebook exercise 3).

Interactive figure: loading…

Sufficiency. A statistic $T(X)$ is sufficient for $\theta$ if the conditional distribution of the data given $T(X)$ does not depend on $\theta$. Once you know $T$, the rest of the data tells you nothing more about $\theta$.

Theorem (factorisation, Fisher–Neyman). $T$ is sufficient for $\theta$ if and only if $p(x \mid \theta) = g(T(x), \theta)\, h(x)$ for some functions $g$ and $h$.

Proof (discrete case)

1. (⇐) Fix $x$ and let $t = T(x)$. Then $P_\theta(X = x \mid T = t) = \dfrac{g(t,\theta)\,h(x)}{\sum_{x' : T(x') = t} g(t,\theta)\,h(x')} = \dfrac{h(x)}{\sum_{x' : T(x') = t} h(x')}$, which does not depend on $\theta$.
2. (⇒) Take $g(t, \theta) = P_\theta(T = t)$ and $h(x) = P(X = x \mid T = T(x))$, which is free of $\theta$ by sufficiency. Then $p(x \mid \theta) = P_\theta(T = T(x))\, P(X = x \mid T = T(x)) = g(T(x), \theta)\, h(x)$. $\blacksquare$

For i.i.d. data from an exponential family, $p(x_1, \dots, x_n \mid \eta) = \big(\prod_i h(x_i)\big)\exp\big(\eta^\top \sum_i T(x_i) - nA(\eta)\big)$, so $\sum_i T(x_i)$ (with $n$) is sufficient. Its size $d$ does not grow with $n$. Coin: the number of heads. Gaussian: $(\sum_i x_i, \sum_i x_i^2)$. Gamma: $(\sum_i \log x_i, \sum_i x_i)$. The same structure gives every exponential family a conjugate prior (lesson 10).

When maximum likelihood fails #

• Zero counts. A new web page variant has $k = 0$ conversions in $n = 20$ visits. The MLE is $\hat p = 0$: "this page never converts". In the natural parameter, $\hat\eta = \log\frac{0}{1} = -\infty$: the maximum is not attained. The same happens for all-zero Poisson counts and for logistic regression on data that a line separates perfectly (the weights grow without bound). A prior fixes this (lesson 10).
• No unique maximum. If the model is not identifiable, a whole set of $\theta$ values gives the same likelihood.
• Small samples. The MLE can be biased (the variance by $(n-1)/n$, the rate by $n/(n-1)$); the bias shrinks like $1/n$.

An unbounded likelihood: Gaussian mixtures

Fit the mixture $\frac12\mathcal N(\mu_1, \sigma_1^2) + \frac12\mathcal N(\mu_2, \sigma_2^2)$ to data $x_1, \dots, x_n$ with at least two different values. Put $\mu_1 = x_1$ and let $\sigma_1 \to 0$. The density at $x_1$ contains the term $\frac{1}{2\sqrt{2\pi}\sigma_1}$, which goes to infinity. Every other point keeps a density of at least $\frac12\mathcal N(x_i; \mu_2, \sigma_2^2) > 0$ from the second component. So $L \to \infty$: the "MLE" is a spike on one data point, a useless model. In practice, EM for Gaussian mixtures needs a floor on the variances or a prior on them (this is MAP estimation, lesson 10). The KL view does not apply here, because no $\theta$ achieves the supremum.

Answers to section 1. (1) Too small, by the factor $(n-1)/n = 0.9$: $\mathbb{E}[\hat\sigma^2] = 0.9\,\sigma^2$. The $1/n$ variance of 10 seeds is 10% low on average. Divide by $n - 1 = 9$ for an unbiased estimate. (2) By invariance, the MLE of $\theta = 1/\lambda$ is $1/\hat\lambda = \bar x$. But $\hat\lambda = 1/\bar x$ is biased upward: $\mathbb{E}[\hat\lambda] = \frac{n}{n-1}\lambda$. The MLE keeps its form under a change of parameter; unbiasedness does not survive a nonlinear map.

4. Worked example #

(a) A kitchen scale. You weigh a 1 kg reference weight four times. The readings minus 1000 g are $x = (2.1,\ -0.3,\ 3.0,\ 0.4)$ g. Model: $x_i \sim \mathcal N(\mu, \sigma^2)$ i.i.d., where $\mu$ is the scale's offset and $\sigma$ its noise.

1. Sufficient statistics. $\sum_i x_i = 2.1 - 0.3 + 3.0 + 0.4 = 5.2$ and $\sum_i x_i^2 = 4.41 + 0.09 + 9.00 + 0.16 = 13.66$. For this model, nothing else in the data matters.
2. MLE of the offset. $\hat\mu = 5.2/4 = 1.3$ g.
3. MLE of the variance by moment matching on $x^2$. $\hat\mu^2 + \hat\sigma^2 = 13.66/4 = 3.415$, so $\hat\sigma^2 = 3.415 - 1.69 = 1.725$ g².
4. Check with the deviations. $x_i - \bar x = (0.8,\ -1.6,\ 1.7,\ -0.9)$. Their squares are $0.64, 2.56, 2.89, 0.81$, with sum $6.9$. And $6.9/4 = 1.725$ ✓.
5. Bias. $\mathbb{E}[\hat\sigma^2] = \frac{3}{4}\sigma^2$: with only 4 readings, the MLE is 25% low on average. The unbiased estimate is $S^2 = 6.9/3 = 2.3$ g².
6. Invariance. The MLE of the noise standard deviation is $\hat\sigma = \sqrt{1.725} \approx 1.313$ g.

(b) Waiting times. A help desk records the times between 5 calls: their sum is $S = 10$ minutes. Model: i.i.d. Exponential($\lambda$).

1. $\ell(\lambda) = n\log\lambda - \lambda S = 5\log\lambda - 10\lambda$.
2. $\ell'(\lambda) = 5/\lambda - 10 = 0$ gives $\hat\lambda = 0.5$ calls per minute. $\ell''(\lambda) = -5/\lambda^2 < 0$: a maximum.
3. Invariance: the MLE of the mean waiting time is $1/\hat\lambda = 2$ minutes $= \bar x$.
4. Bias: $\mathbb{E}[\hat\lambda] = \frac{n}{n-1}\lambda = 1.25\,\lambda$. With 5 observations, the rate estimate is 25% too high on average. $\frac{n-1}{n}\hat\lambda = \frac{n-1}{S} = 0.4$ is unbiased.

(c) The coin in natural form. For $k = 7$ heads in $n = 10$ flips, the log-likelihood in the log-odds $\eta$ is $\ell(\eta) = 7\eta - 10\log(1 + e^\eta)$.

1. $\ell'(\eta) = 7 - 10\,\sigma(\eta) = 0$ gives $\sigma(\hat\eta) = 0.7$: the model mean matches the sample mean.
2. $\hat\eta = \log\frac{0.7}{0.3} \approx 0.847$, which is $\log\frac{\hat p}{1-\hat p}$ (invariance).
3. $\ell''(\eta) = -10\,\sigma(\eta)(1 - \sigma(\eta)) = -2.1$ at $\hat\eta$, negative for every $\eta$: $\ell$ is concave in $\eta$.

5. In code #

NumPy and SciPy check the worked example, the bias of the variance estimators and the gamma moment matching. The notebook has these cells and the exercises.

import numpy as np
from scipy import stats, special

# (a) kitchen scale: MLE mean, MLE variance (ddof=0) and the unbiased S^2 (ddof=1)
x = np.array([2.1, -0.3, 3.0, 0.4])
print(x.mean(), x.var(), round(x.var(ddof=1), 6))      # 1.3 1.725 2.3

# predict-first 1: 10 seeds; the 1/n estimate averages 0.9 * sigma^2
rng = np.random.default_rng(1)
runs = rng.normal(0.80, 0.02, size=(200_000, 10))      # accuracy 0.80, sd 0.02
print(round(runs.var(axis=1).mean() / 0.02**2, 3))    # 0.898 (exact: 0.9)

# bias of 1/xbar for exponential data, n = 10, lambda = 1 (exact: mean 1.111, MSE 0.167)
lam_hat = 1 / rng.exponential(1.0, size=(400_000, 10)).mean(axis=1)
print(round(lam_hat.mean(), 3), round(((lam_hat - 1) ** 2).mean(), 3))   # 1.111 0.166

# gamma MLE = moment matching on T = (log x, x); NumPy and SciPy use scale = 1 / rate
g = np.random.default_rng(3).gamma(shape=2.0, scale=1 / 0.5, size=200)
k_hat, _, scale_hat = stats.gamma.fit(g, floc=0)
b_hat = 1 / scale_hat
print(np.isclose(special.digamma(k_hat) - np.log(b_hat), np.log(g).mean()))   # True
print(np.isclose(k_hat / b_hat, g.mean()))                                     # True

The printed Monte Carlo values move a little with the seed; their errors are about $\pm 0.002$ here. The exact values come from the formulas of section 3.

6. Where it shows up #

• Every supervised loss. Mean squared error, binary and softmax cross-entropy are Gaussian, Bernoulli and categorical NLLs. For the categorical, the logits $z$ are the natural parameters, $A(z) = \log\sum_k e^{z_k}$ and $\nabla A = \operatorname{softmax}(z)$. So the gradient of the loss with respect to the logits is "softmax minus one-hot", the rule of calculus 06.
• Generalised linear models. With the canonical link, $\eta_i = w^\top x_i$, and the NLL is $\sum_i \big[A(w^\top x_i) - y_i\, w^\top x_i\big] + \text{const}$: convex in $w$, with gradient $\sum_i (\hat y_i - y_i)\, x_i$, where $\hat y_i = A'(w^\top x_i)$. Linear, logistic and Poisson regression have the same form.
• Language models. Pretraining minimises the average NLL of the next token: the MLE of a categorical model. Perplexity is $e^{\text{average NLL}}$.
• Models that predict their own uncertainty. Heteroscedastic regression, probabilistic weather models and surrogates of simulators predict a mean and a variance and minimise the Gaussian NLL.
• Science. Gamma fits of rainfall amounts (the Standardized Precipitation Index), exponential fits of the magnitudes of earthquakes above a threshold (Aki's estimator of the Gutenberg–Richter $b$-value is an exponential MLE), Poisson rates of photon counts in astronomy.
• Evaluation. Reporting the variance of a metric over a handful of seeds or folds with the $1/n$ formula underestimates it; with 3 seeds the factor is $2/3$.

7. Common mistakes #

• Treating the likelihood as a probability of $\theta$. $L(\theta)$ is the probability of the data. It does not integrate to 1 over $\theta$, and "the probability that $p = 0.7$" needs a prior (lesson 10).
• Thinking unbiased means best. Compare MSEs. The $1/n$ and $1/(n+1)$ variance estimators beat the unbiased $S^2$ in MSE for Gaussian data.
• Expecting unbiasedness to survive a transformation. $S^2$ is unbiased for $\sigma^2$, but $S$ is biased for $\sigma$ (Jensen). The MLE does survive transformations (invariance).
• Differentiating $A$ in the wrong parameters. $\nabla A = \mathbb{E}[T]$ holds for the natural parameter $\eta$, not for $(\mu, \sigma)$ or $(k, \beta)$.
• Mixing up rate and scale. These lessons write the exponential and gamma with the rate. NumPy and SciPy use the scale $1/\beta$.
• Assuming the MLE always exists. It does not for all-zero counts, separable data in logistic regression, constant data in a Gaussian or gamma fit, or Gaussian mixtures.
• Maximising the likelihood with a product. Multiply 1000 probabilities and you get 0 in floating point. Always sum log-probabilities.

8. Check yourself #

1. Write the exponential distribution in exponential-family form, and get its mean and variance from $A$.
   Answer

   $p(x) = \lambda e^{-\lambda x} = \exp(-\lambda x + \log\lambda)$ on $x \ge 0$. So $\eta = -\lambda$, $T(x) = x$, $h(x) = 1$ on $x \ge 0$, and $A(\eta) = -\log\lambda = -\log(-\eta)$. Then $A'(\eta) = -1/\eta = 1/\lambda = \mathbb{E}[X]$ and $A''(\eta) = 1/\eta^2 = 1/\lambda^2 = \operatorname{Var}(X)$.
2. A website gets $x = (5, 7, 4, 6, 3)$ sign-ups in five hours, modelled as i.i.d. Poisson($\lambda$). Give the MLE, check that it is a maximum, and give the MLE of $P(\text{no sign-up in an hour}) = e^{-\lambda}$.
   Answer

   $\ell(\lambda) = 25\log\lambda - 5\lambda + \text{const}$, so $\ell'(\lambda) = 25/\lambda - 5 = 0$ gives $\hat\lambda = 5$. $\ell''(\lambda) = -25/\lambda^2 < 0$: a maximum. By invariance, the MLE of $e^{-\lambda}$ is $e^{-5} \approx 0.0067$.
3. $S^2$ is unbiased for $\sigma^2$. Is $S$ unbiased for $\sigma$?
   Answer

   No. The square root is strictly concave, so by Jensen $\mathbb{E}[S] < \sqrt{\mathbb{E}[S^2]} = \sigma$ (strict, because $S$ is not constant). For Gaussian data with $n = 5$, $\mathbb{E}[S] \approx 0.94\,\sigma$.
4. Which of these estimators of a mean are consistent? (a) $X_1$; (b) $\bar X + 5/n$; (c) $\bar X$ for data from a Cauchy distribution centred at $m$, as an estimator of $m$.
   Answer

   (a) No: it is unbiased, but its distribution does not change with $n$. (b) Yes: the bias $5/n \to 0$ and the variance $\sigma^2/n \to 0$, so the MSE goes to 0. (c) No: the mean of $n$ i.i.d. Cauchy variables has the same Cauchy distribution for every $n$. The Cauchy has no mean, so the law of large numbers does not apply. The sample median is consistent for $m$.
5. Show that the binary cross-entropy of a logistic model, $-\sum_i [y_i\log\sigma(z_i) + (1 - y_i)\log(1 - \sigma(z_i))]$, is the negative Bernoulli log-likelihood, and give its derivative with respect to one logit $z_i$.
   Answer

   The Bernoulli PMF with $p_i = \sigma(z_i)$ is $p_i^{y_i}(1 - p_i)^{1 - y_i}$; its log is $y_i\log p_i + (1 - y_i)\log(1 - p_i)$. Sum and negate: that is the loss. In natural form, the term is $-(y_i z_i - \log(1 + e^{z_i}))$, so the derivative is $\sigma(z_i) - y_i$: prediction minus target.
6. What are the sufficient statistics of a gamma sample? Why does a gamma fit fail for constant data?
   Answer

   $T(x) = (\log x, x)$, so $(\sum_i \log x_i, \sum_i x_i)$, or the geometric and the arithmetic mean. The shape equation is $\log k - \psi(k) = \log\bar x - \overline{\log x}$. For constant data the right side is 0 (the two means are equal), and $\log k - \psi(k) > 0$ for every $k > 0$, decreasing to 0 only as $k \to \infty$. So no MLE exists: the fit wants a spike of infinite shape.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the likelihood $L(\theta)$?

   Show answer

2. MLE of a Bernoulli $p$ from $k$ successes in $n$ trials?

   Show answer

3. Gaussian MLEs of $\mu$ and $\sigma^2$?

   Show answer

4. State the bias–variance decomposition of the MSE.

   Show answer

5. What is $\mathbb{E}[\hat\sigma^2_{\text{MLE}}]$ for $n$ Gaussian observations?

   Show answer

6. What does consistency of an estimator mean, and what is a sufficient condition?

   Show answer

7. Why does maximum likelihood minimise a KL divergence?

   Show answer

8. Write the general exponential-family form.

   Show answer

9. What are $\nabla A(\eta)$ and $\nabla^2 A(\eta)$ in an exponential family?

   Show answer

10. What is the invariance property of the MLE?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Derive the MLE of a Gaussian. Why does the unbiased variance divide by $n - 1$?" — $\hat\mu = \bar x$ and $\hat\sigma^2 = \frac1n\sum_i (x_i - \bar x)^2$; because $\mathbb{E}\sum_i (X_i - \bar X)^2 = (n - 1)\sigma^2$: the deviations are measured from $\bar X$, which is fitted to the same data.
• "Mean squared error and cross-entropy: which likelihoods?" — Gaussian with a fixed variance; Bernoulli or categorical. Minimising them is maximum likelihood, which minimises $\mathrm{KL}(p_{\text{data}} \,\|\, p_\theta)$.
• "Is an unbiased estimator always better?" — No. $\text{MSE} = \text{bias}^2 + \text{variance}$; shrinkage (ridge, MAP, James–Stein) trades bias for variance.
• "Why is the logistic-regression loss convex?" — It is an exponential-family NLL with a natural parameter linear in $w$: the convex $A(w^\top x_i)$ minus a linear term.
• "Your logistic regression weights blow up. Why?" — The data are linearly separable, so the likelihood keeps increasing as $\|w\| \to \infty$; the MLE does not exist. Add L2 regularisation (a prior).
• "What is a sufficient statistic? Give one for a Gaussian." — A summary that carries all the information about $\theta$; $(\sum_i x_i, \sum_i x_i^2)$.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. The MLE of a Poisson rate by a grid search over the log-likelihood, compared with the closed form $\bar x$.
2. Bias, variance and MSE by simulation (mc_bias_variance), applied to the Gaussian variance estimators with $1/(n-1)$, $1/n$ and $1/(n+1)$.
3. The gamma MLE by Newton's method on $\log k - \psi(k) = \log\bar x - \overline{\log x}$, compared with scipy.stats.gamma.fit(x, floc=0).
4. Bias, variance and MSE of $1/\bar x$ for exponential data with $n = 10$, against the closed form.
5. The cross-entropy loss of a logistic model equals the Bernoulli NLL; check it and its gradient numerically.

11. Watch and read #

Maximum Likelihood, clearly explained!!!StatQuest with Josh Starmer · 6 min

The idea of the likelihood as a curve over the parameter, with a fitted distribution moving over the data.

Maximum Likelihood For the Normal Distribution, step-by-step!!!StatQuest with Josh Starmer · 20 min

The Gaussian MLE derived slowly, with the log-likelihood and both partial derivatives.

5. Maximum Likelihood Estimation (cont.)MIT OpenCourseWare · Philippe Rigollet, 18.650 Statistics for Applications · 77 min

The MLE as an estimate of the KL divergence, with Bernoulli, Poisson and Gaussian examples.

• Kevin Murphy, Probabilistic Machine Learning: An Introduction, free PDF: §4.2 (maximum likelihood; §4.2.2 is the KL justification) and §4.7.6 (the bias–variance trade-off).
• Kevin Murphy, Probabilistic Machine Learning: Advanced Topics, free PDF: §2.4 (the exponential family: log-partition function, moments, MLE).
• Philippe Rigollet, MIT 18.650 Statistics for Applications, lecture slides and videos (lectures 4–6: maximum likelihood and the method of moments).
• Thomas Minka, Estimating a Gamma distribution (2002): the gamma MLE, a fast Newton update and a starting value (it uses the scale $1/\beta$).