Maximum Likelihood Estimation

Let $\{f(x; \theta) : \theta \in \Theta\}$ be a parametric family of densities (continuous case) or PMFs (discrete case), and let $\mathbf{x} = (x_1, \ldots, x_n)$ be an observed sample. The likelihood function is the joint density or PMF evaluated at the observed sample, viewed as a function of the parameter $\theta$ with the data held fixed:

When the data are understood we write $L(\theta)$ for short.

The log-likelihood is

It is the sum (not the product) of per-observation log-densities, and it has the same maximizer as $L(\theta)$ because $\log$ is strictly increasing.

The switch from $L$ to $\ell$ is not cosmetic. Three things happen simultaneously, and each is load-bearing.

(1) Products become sums. For iid data $L(\theta) = \prod_i f(X_i; \theta)$ is a product of $n$ terms; $\ell(\theta) = \sum_i \log f(X_i; \theta)$ is a sum. Sums are what every limit theorem in probability speaks to — the law of large numbers for sample averages, the central limit theorem for standardized sums. The transition $L \to \ell$ is what lets consistency and asymptotic normality inherit directly from the LLN and CLT applied to $\log f(X; \theta)$.

(2) Concavity (often). For exponential families in the natural parameter, $\ell(\theta)$ is strictly concave: it has a unique maximizer, and any local maximum is the global one. The likelihood $L(\theta) = e^{\ell(\theta)}$ is not concave, just log-concave — a strictly weaker property. Working in log-space is what makes uniqueness and Newton-Raphson convergence arguments tractable.

(3) Numerical stability. Products of many small probabilities underflow to zero in floating point. A sample of $n = 1000$ Normal observations produces an $L$ that is effectively zero — but $\ell$ is a perfectly ordinary number on the order of $-n$. Anyone who has fit a model to more than a few hundred data points has implicitly used $\ell$ for this reason alone.

Let $X_1, \ldots, X_n$ be iid $\mathcal{N}(\mu, \sigma^2)$ with $\sigma^2$ known. The per-observation density is $f(x; \mu) = (2\pi\sigma^2)^{-1/2} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$. Taking the product over $i$ and then the logarithm:

The first term does not depend on $\mu$, so we can drop it when maximizing; what remains is the negative sum of squared residuals. Maximizing $\ell(\mu)$ is therefore the same as minimizing $\sum_i (X_i - \mu)^2$. Least squares and maximum likelihood, for Normal-error models, are the same procedure seen from two angles — one statistical, one geometric.

Let $X_1, \ldots, X_n$ be iid $\operatorname{Bernoulli}(p)$. The PMF is $f(x; p) = p^x (1-p)^{1-x}$ for $x \in \{0, 1\}$. Writing $k = \sum_i X_i$ for the number of successes,

The likelihood is a polynomial in $p$ of degree $n$; the log-likelihood is linear in $\log p$ and $\log(1-p)$. The shape of $\ell(p)$ — a single peak somewhere between $0$ and $1$, concave — is already visible in the formula without any calculus.

Let $X_1, \ldots, X_n$ be iid $\operatorname{Exponential}(\lambda)$ with density $f(x; \lambda) = \lambda e^{-\lambda x}$ for $x > 0$. Then

Again a clean closed-form log-likelihood, linear in $\lambda$ apart from the $n \log \lambda$ growth term. For any sample with $\bar{X}_n > 0$, $\ell$ is strictly concave in $\lambda$ (its second derivative is $-n/\lambda^2 < 0$), so the maximizer is unique.

The maximum-likelihood estimator of $\theta$ is any value $\hat{\theta}_n \in \Theta$ that maximizes $L(\theta)$ — equivalently, that maximizes $\ell(\theta)$:

When the maximizer is unique, we write $\hat{\theta}_n = \arg\max_\theta \ell(\theta)$. When multiple maxima exist — e.g., on the boundary of $\Theta$, or at non-smooth points — we pick any of them; the finite-sample behavior can differ, but the asymptotic theory is unaffected.

Suppose $\Theta \subseteq \mathbb{R}$ is an open interval, $\ell(\theta)$ is differentiable on $\Theta$, and the global maximizer $\hat{\theta}_n$ lies in the interior of $\Theta$. Then $\hat{\theta}_n$ satisfies the score equation

where $s(\theta; x) = \partial_\theta \log f(x; \theta)$ is the per-observation score.

From Example 1, $\ell(\mu) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_i (X_i - \mu)^2$. Differentiating,

Setting the derivative to zero and dividing by $n$ gives $\sum_i X_i = n\mu$, so

The second derivative $-n/\sigma^2 < 0$ confirms this is a maximum, not a minimum. The MLE of the Normal mean is the sample mean. Every other Topic 13 property of $\bar{X}_n$ — unbiasedness, variance $\sigma^2/n$, asymptotic normality with rate $\sqrt{n}$, efficiency — is now a property of the MLE.

From Example 2, $\ell(p) = k \log p + (n-k) \log(1-p)$ with $k = \sum_i X_i$. Differentiating,

Setting this to zero: $k(1-p) = (n-k)p$, so $k - kp = np - kp$, and therefore $k = np$, giving

The sample proportion is the MLE — unsurprising in hindsight, and intuitive. The second derivative is $-k/p^2 - (n-k)/(1-p)^2 < 0$ for $p \in (0, 1)$, confirming a maximum. At the boundary $k = 0$ or $k = n$ the MLE is $0$ or $1$, at the edge of the parameter space — a clean example of a boundary solution.

From Example 3, $\ell(\lambda) = n \log \lambda - \lambda n \bar{X}_n$. Differentiating,

Setting this to zero gives $\lambda = 1/\bar{X}_n$:

Note that the MLE is a nonlinear function of $\bar{X}_n$. Because $1/x$ is convex, Jensen’s inequality implies $\mathbb{E}[1/\bar{X}_n] > 1/\mathbb{E}[\bar{X}_n] = \lambda$, so the MLE is biased for finite $n$. We will recover bias-free behavior asymptotically, but this is a useful reminder that MLEs are not generally unbiased; they are consistent and asymptotically efficient, not unbiased in finite samples.

For iid $\operatorname{Poisson}(\lambda)$ data with PMF $f(x; \lambda) = \lambda^x e^{-\lambda}/x!$,

The final $\log(X_i!)$ term does not depend on $\lambda$. Differentiating the rest,

and setting this to zero gives $\hat{\lambda}_{\text{MLE}} = \bar{X}_n$. The Poisson MLE is the sample mean — the same estimator as for the Normal mean and the Bernoulli proportion. This is not a coincidence: in each case the natural sufficient statistic is $\sum_i X_i$, and the MLE is its normalized version. Exponential families make this pattern universal (Exponential Families, §7.8).

The four examples above all give the MLE as a clean function of the sample mean. That is a feature of their exponential-family structure, not a general property of MLE. For most parametric families — Gamma shape, beta shape, Weibull, Cauchy location, logistic regression, hidden Markov models, neural networks — the score equation has no closed-form solution, and the MLE must be computed numerically. Section 14.7 develops Newton-Raphson for exactly this situation: write down the score and the observed Fisher information, pick a starting value, iterate $\theta^{(t+1)} = \theta^{(t)} + S(\theta^{(t)}) / J(\theta^{(t)})$, and stop when the steps become small. Every optimization method in modern machine learning — SGD, Adam, L-BFGS — is a variation on this theme applied to the log-likelihood or a regularized variant.

Let $\hat{\theta}_n$ be the MLE of $\theta$, and let $g : \Theta \to \mathcal{G}$ be any function (not necessarily one-to-one). Then the MLE of $\phi = g(\theta)$ is

If $\hat{\sigma}^2_{\text{MLE}} = \frac{1}{n}\sum_i (X_i - \bar{X}_n)^2$ is the MLE of the Normal variance $\sigma^2$ (as derived from the two-parameter Normal log-likelihood by zeroing the score in $\sigma^2$), then by invariance the MLE of the standard deviation $\sigma = \sqrt{\sigma^2}$ is

No separate maximization is needed. Notice that $\hat{\sigma}_{\text{MLE}}$ is biased even though $\hat{\sigma}^2_{\text{MLE}}$ would be unbiased if we used $n-1$ instead of $n$; invariance preserves the MLE but does not preserve unbiasedness under nonlinear transformations. The two properties simply track different things.

In Bernoulli data, the MLE of $p$ is $\hat{p} = \bar{X}_n$. The odds are $\phi(p) = p/(1-p)$, a strictly increasing transformation on $(0, 1)$ (so the one-to-one case applies). By invariance, the MLE of the odds is

The odds are the natural parameter of the Bernoulli family: $\eta = \log\!\frac{p}{1-p}$ is the logit link that logistic regression uses. Invariance is what makes “logistic regression on $\eta$” equivalent to “MLE of $p$ then take the log-odds” — the two routes agree at the estimate.

Invariance transports the point estimate under $g$, but it does not transport the sampling distribution. If $\sqrt{n}(\hat{\theta}_n - \theta_0) \overset{d}{\to} \mathcal{N}(0, 1/I(\theta_0))$, what happens to $\sqrt{n}(g(\hat{\theta}_n) - g(\theta_0))$? Theorem 7 in Topic 13 (delta method) gives the answer:

So the asymptotic variance scales by $[g'(\theta_0)]^2$. A consequence: the Fisher information transforms tensorially, and efficient estimators in one parameterization remain efficient under smooth reparameterization. Invariance is about the point; the delta method is about the spread.

The family $\{f(x; \theta) : \theta \in \Theta\}$ is Kullback–Leibler identifiable at $\theta_0$ if for every $\theta \neq \theta_0$ in $\Theta$,

Equivalently: $\theta = \theta_0$ is the unique maximizer of $\theta \mapsto \mathbb{E}_{\theta_0}[\log f(X; \theta)]$.

Let $X_1, X_2, \ldots$ be iid from $f(x; \theta_0)$ where $\theta_0 \in \Theta \subseteq \mathbb{R}$. Assume (i) $\Theta$ is compact, (ii) the family is KL-identifiable at $\theta_0$, and (iii) $\theta \mapsto \log f(x; \theta)$ is continuous in $\theta$ for every $x$ and dominated by an integrable function. Then the MLE $\hat{\theta}_n$ is weakly consistent:

The three conditions in Theorem 3 — compact $\Theta$, KL-identifiability, continuity with domination — are called regularity conditions. They are standard in most parametric problems, and they can be relaxed in a number of directions:

• Compactness of $\Theta$: Usually replaced by “the MLE eventually lies in a compact subset of $\Theta$ with probability one,” which holds under consistency of a coarser preliminary estimator. For open $\Theta$ (e.g., $\Theta = \mathbb{R}$ for the Normal mean), this is almost always true, but the proof becomes messier.
• Identifiability: Non-negotiable. Label-switching in mixture models, equivalent parameterizations in neural networks, and unobservable components in latent-variable models all violate identifiability in different ways. The MLE in these settings is consistent only for the identified quotient — e.g., the set of mixture distributions, not the parameter vector of component labels.
• Continuity and domination: Continuity of $\log f$ in $\theta$ is almost always immediate. Domination by an integrable function (so expectations are finite and uniform LLN applies) can fail for heavy-tailed families, which is why the MLE story becomes more delicate for Cauchy and similar distributions.

When the conditions hold, they buy you consistency for free — no explicit convergence rate, no finite-sample guarantees, but the estimator does eventually home in on the truth. Rates and finite-sample behavior are the job of asymptotic normality (§14.5) and large-deviations-type bounds (Topic 12).

Let $X_1, X_2, \ldots$ be iid from $f(x; \theta_0)$ on an open interval $\Theta$, with $\hat\theta_n$ the MLE. Assume the conditions of Theorem 3 for consistency, plus (iv) $\theta \mapsto \log f(x; \theta)$ is twice continuously differentiable in a neighborhood of $\theta_0$ with dominated second derivative, and (v) the Fisher information $I(\theta_0) \in (0, \infty)$. Then

For Normal$(\mu, \sigma^2)$ with $\sigma^2$ known and $\hat\mu_n = \bar{X}_n$, $I(\mu) = 1/\sigma^2$ (Example 10 of Topic 13). Theorem 4 gives

This is the classical CLT for the sample mean, now recovered as a special case of the general MLE asymptotic-normality theorem. For Normal data, the “asymptotic” normality is actually exact at every finite $n$: $\bar{X}_n \sim \mathcal{N}(\mu, \sigma^2/n)$ for any $n \geq 1$, not just in the limit. The MLE framework gives one unified result; for this particular family the result happens to hold non-asymptotically.

For Bernoulli$(p)$ with $\hat p_n = \bar{X}_n$, the Fisher information is $I(p) = 1/(p(1-p))$, so $1/I(p) = p(1-p)$. Theorem 4 gives

This is the two-sided de Moivre–Laplace theorem in disguise. For a 95% Wald confidence interval we replace $p$ by $\hat p$ in the asymptotic variance:

This is the formula every introductory statistics class derives ad hoc — it is the MLE Wald interval applied to the Bernoulli family. The interval has known coverage issues near $p = 0$ or $p = 1$; the Wilson score interval and the Agresti–Coull interval are alternatives. What the MLE machinery gives us is the default construction; better intervals require more careful analysis of finite-sample behavior.

The arrow $\overset{d}{\to}$ is deceptively clean. In finite samples, three things can go wrong:

(1) Slow convergence near boundaries. For Bernoulli with true $p$ near $0$ or $1$, the sampling distribution of $\hat p$ is very skewed and discrete; the Wald interval undercovers dramatically. The convergence rate is still $\sqrt n$, but the hidden constants are huge near boundary points.

(2) Bias at second order. The MLE is consistent but generally not unbiased; the bias shrinks at rate $1/n$, so it is negligible asymptotically but can matter at moderate sample sizes. The classical “Bartlett correction” adjusts likelihood-ratio tests for this second-order bias.

(3) Asymptotic efficiency ≠ finite-sample optimality. At fixed $n$, a biased estimator with smaller variance (e.g., James–Stein in $d \geq 3$, or a shrinkage estimator) can dominate the MLE in MSE. Asymptotic efficiency says the MLE is optimal as $n \to \infty$; for small $n$ you may well want to look elsewhere. Topic 13 §13.8 developed this point carefully.

A healthy reading of Theorem 4 is: “for large enough $n$, with well-behaved parameters, you will not be led astray by the Gaussian approximation.” The phrases “large enough” and “well-behaved” are load-bearing and should be checked by simulation whenever $n$ is less than a few hundred or the parameter is near a boundary.

Under the conditions of Theorem 4, the MLE achieves the Cramér–Rao asymptotic variance: its asymptotic variance $1/I(\theta_0)$ equals the Cramér–Rao lower bound for the class of regular estimators. In particular, for any other regular asymptotically-normal estimator $\tilde\theta_n$ with $\sqrt{n}(\tilde\theta_n - \theta_0) \overset{d}{\to} \mathcal{N}(0, v(\theta_0))$,

Line up the three workhorse examples:

All three MLEs achieve the CRLB in the limit; the Normal mean achieves it exactly at every finite $n$. Efficiency is not specific to exponential families — it is specific to the MLE — but exponential families make the arithmetic particularly clean because the score is linear in the sufficient statistic.

Theorem 5 says no regular estimator can beat the MLE asymptotically. The catch is the word regular. Joseph Hodges in 1951 constructed the following superefficient estimator: shrink $\bar{X}_n$ all the way to $0$ whenever $|\bar{X}_n| < n^{-1/4}$, and use $\bar{X}_n$ otherwise. At $\theta_0 = 0$ the estimator has zero asymptotic variance; at every other $\theta_0$ it behaves like $\bar{X}_n$.

This does not violate Theorem 5 because the Hodges estimator is irregular at $\theta_0 = 0$ — the sampling distribution does not converge uniformly in $\theta_0$ near zero. More fundamentally, the maximum risk over a shrinking neighborhood of $\theta_0 = 0$ diverges: the estimator pays for its superefficiency at one point with enormous variance at nearby points. Le Cam’s convolution theorem makes this precise: any estimator that is asymptotically better than the MLE at one point is asymptotically worse at nearby points, in a specific technical sense.

For practical purposes, the lesson is that superefficiency is a local, fragile phenomenon. The MLE is the default for good reason, and the escape routes — shrinkage, empirical Bayes, Stein-type estimators — deliver their gains through bias, regularization, or prior information, not through asymptotic efficiency.

Given an iterate $\theta^{(t)}$, the Newton-Raphson update linearizes the score equation around $\theta^{(t)}$ and solves the linear equation:

where $J(\theta) = -\ell_n''(\theta) = -S_n'(\theta)$ is the observed Fisher information. The update is applied until $|\theta^{(t+1)} - \theta^{(t)}|$ falls below a tolerance.

Suppose $\ell_n$ is three times continuously differentiable in a neighborhood of the MLE $\hat\theta$, and $J(\hat\theta) > 0$. Then for any starting value $\theta^{(0)}$ sufficiently close to $\hat\theta$, the Newton-Raphson iterates satisfy

for some constant $C > 0$ depending on $\hat\theta$ and the third derivative of $\ell_n$. Equivalently, the number of correct digits roughly doubles per iteration.

For iid Gamma$(\alpha, \beta)$ data with $\beta$ known, the log-likelihood is

Differentiate once to get the score and once more to get the Hessian:

where $\psi(\alpha) = \Gamma'(\alpha)/\Gamma(\alpha)$ is the digamma function and $\psi'(\alpha)$ is the trigamma function. The observed information is $J(\alpha) = n\psi'(\alpha) > 0$, so Newton-Raphson is

A standard initial guess is the method-of-moments estimate $\alpha^{(0)} = \beta \bar{X}_n$, since $\mathbb{E}[X] = \alpha/\beta$. From this start, convergence typically takes $4$–$6$ iterations to machine precision. No closed-form solution exists because $\psi^{-1}$ has no closed form — but we do not need one; Newton-Raphson handles it numerically.