The Method of Maximum Likelihood

Derivation

The Method of Maximum Likelihood

D-119 Home PU-106 Threads chance Depends on Expectation, Variance and Their Algebra

Statement

Given a dataset x = (x₁, …, x_n) modelled as independent draws from a density (or mass function) f(x; θ) indexed by a parameter θ, the likelihood is the joint density read as a function of θ with the data held fixed, L(θ) = ∏_i f(x_i; θ). The maximum-likelihood estimator θ̂ maximises L, equivalently the log-likelihood ℓ(θ) = ∑_i ln f(x_i; θ), and for a regular interior maximum it solves the estimating (score) equation ∑_i ∂ ln f(x_i; θ)/∂θ = 0.

Why it matters

Maximum likelihood is the default engine of statistical inference in physics: it turns a probabilistic model plus a finite dataset into a single best-fit parameter, with a principled recipe rather than an ad-hoc curve fit. Least-squares fitting, the sample mean, and the sample variance all fall out of it as special cases.

It also comes with a companion error theory. The curvature of ℓ at its peak — the observed information — gives the estimator's variance, so maximum likelihood delivers not just θ̂ but its uncertainty, and under regularity it attains the smallest variance any unbiased estimator can have (the Cramér–Rao bound).

Assumptions

The data are independent draws from the model.If the observations are correlated, the joint density does not factorise into ∏_i f(x_i; θ) and the log-likelihood is not a simple sum; the score equation below acquires cross terms.

The model family contains a well-defined f(x; θ) that is a proper normalised density for every θ.If f is unnormalised or the parametrisation is not identifiable (two different θ giving the same distribution), the maximiser is not unique or not meaningful.

The likelihood is a smooth function of θ on an open parameter set (regularity).If ℓ is non-differentiable or the true value sits on the boundary of the allowed region, setting the derivative to zero can miss the maximum entirely.

The support of f does not depend on θ.If the range of allowed x shifts with θ (as in uniform(0, θ)), differentiation under the integral sign fails, the score need not have zero mean, and the stationary-point argument collapses — the MLE lives at a boundary instead.

Differentiation and integration (or summation) may be interchanged.If ∂/∂θ cannot be moved inside ∫f dx = 1, the key identity E[score] = 0 that underlies consistency and the information calculation no longer holds.

Derivation

p(x₁, …, x_n; θ) = ∏_i=1ⁿ f(x_i; θ)

Independence factorises the joint density into a product of marginals. A

L(θ) ≡ p(x; θ) = ∏_i=1ⁿ f(x_i; θ)

Reinterpret the same expression: the data x are now fixed constants and θ is the variable. This relabelling defines the likelihood. A

ℓ(θ) ≡ ln L(θ) = ∑_i=1ⁿ ln f(x_i; θ)

Take the logarithm. Because ln is strictly increasing, arg max_θ L = arg max_θ ℓ, so the maximiser is unchanged; the product becomes a sum, which is easier to differentiate and numerically stable. B

at a maximum: dℓ/dθ |_θ=θ̂ = 0

A smooth function attains an interior maximum only where its first derivative vanishes (Fermat's stationarity condition). A

dℓ/dθ = ∑_i=1ⁿ ∂/∂θ ln f(x_i; θ)

Differentiate the sum term by term (linearity of d/dθ). A

define the score s(x; θ) ≡ ∂ ln f/∂θ = (∂f/∂θ) / f

Name the per-observation contribution; the second form uses the chain rule d(ln f) = df/f. B

∑_i=1ⁿ s(x_i; θ̂) = 0

Combine steps 4–6: the estimating equation. Solve for θ̂. B

require d²ℓ/dθ² |_θ̂ < 0

The second-order condition distinguishes a maximum from a minimum or inflection. Its negative, −d²ℓ/dθ², is the observed information J(θ̂), whose reciprocal estimates var(θ̂). C

(vector case) ∇_θ ℓ = 0, i.e. ∑_i ∂ ln f(x_i; θ)/∂θ_k = 0 for each k

For a parameter vector θ = (θ₁, …, θ_p) the estimating equation becomes one scalar equation per component; the maximum requires a negative-definite Hessian. C

Result

∑_i=1ⁿ ∂/∂θ ln f(x_i; θ) = 0 at θ = θ̂

Reading. The maximum-likelihood estimate is the parameter value at which the total score — the summed sensitivity of each data point's log-density to θ — cancels to zero. Each observation "pulls" θ in the direction that would make that observation more probable; the MLE is the balance point where all the pulls sum to nothing. A negative second derivative confirms it is a peak, and its magnitude sets the estimator's precision.

Units check. ln f is dimensionless (or carries only the fixed [1/x] of a density, independent of θ), so ∂ ln f/∂θ has units of [θ]⁻¹. The sum of n such terms still has units [θ]⁻¹, and setting a quantity of units [θ]⁻¹ to zero is dimensionally consistent. The observed information d²ℓ/dθ² then has units [θ]⁻², so its reciprocal correctly has units [θ]² — a variance.

Limiting cases

Large n: by the law of large numbers ℓ(θ)/n → E[ln f(X; θ)], whose maximiser is the true θ₀ (minimum KL divergence). Hence θ̂ → θ₀ (consistency).
Gaussian model with known variance: the score equation reduces to ∑(x_i − μ) = 0, giving μ̂ = x̄, the ordinary sample mean.
Exponential-family models: the score is linear in the sufficient statistic, so the estimating equation reads "model mean = sample mean" and MLE coincides with the method of moments.
Flat prior: the Bayesian posterior mode (MAP estimate) equals the MLE, since posterior ∝ likelihood when the prior is constant.
Independent-Gaussian errors with common variance: maximising ℓ is identical to minimising ∑(x_i − model)² — least squares emerges as a special case.

Breaks when

The likelihood is unbounded. In a Gaussian mixture, letting one component's variance σ² → 0 while its mean sits on a data point sends L → ∞. No global maximum exists; the score equation has spurious solutions and naive optimisation diverges.
The support depends on the parameter. For uniform(0, θ), L(θ) = θ⁻ⁿ for θ ≥ max x_i and 0 below it. It is decreasing wherever it is nonzero, so dℓ/dθ = −n/θ is never zero; the MLE is the boundary value θ̂ = max x_i, not a stationary point.
The true parameter is on the boundary. Testing whether a variance component is zero, or a mixing weight is zero, places θ₀ on the edge of the allowed set; the derivative need not vanish there and the usual asymptotic distribution of θ̂ is no longer Gaussian.
The sample is too small. The MLE can be badly biased (e.g. the MLE of Gaussian variance divides by n, not n−1, underestimating σ²) and the curvature-based error bars, which are asymptotic, are unreliable.

Failure modes

Maximising L directly and underflowing. A product of thousands of small densities rounds to zero in floating point; always work with the log-likelihood sum.
Skipping the second-order check. Solving dℓ/dθ = 0 and reporting the root without verifying d²ℓ/dθ² < 0 — the stationary point may be a minimum or a saddle.
Confusing likelihood with probability. Writing "the probability that θ equals …" — L(θ) is not a density in θ, does not integrate to 1 over θ, and cannot be interpreted as P(θ) without a prior.
Applying the score equation to a θ-dependent support. Differentiating uniform(0, θ) and setting it to zero, missing that the answer is the boundary max x_i.
Reporting the MLE variance as unbiased. Using σ̂² = (1/n)∑(x_i−x̄)² where an unbiased estimate needs 1/(n−1).
Trusting a single local optimum. Multimodal likelihoods (mixtures, latent-variable models) need multiple starting points; a single gradient ascent may land on the wrong peak.

Discussion

The conceptual pivot of the method is the reinterpretation in step 2. The very same expression ∏ f(x_i; θ) is a probability density when θ is fixed and x varies (it integrates to 1 over the data space), but a likelihood when x is fixed and θ varies (it need not integrate to anything in particular over θ). Fisher's insight was that this second reading is the natural quantity to optimise: we choose the parameter under which the data we actually observed were most probable.

The score s(X; θ) has a beautiful property. Because ∫f dx = 1 for every θ, differentiating gives ∫∂f/∂θ dx = 0, and rewriting ∂f/∂θ = f ∂ln f/∂θ shows E[s] = 0 at the true θ₀. The MLE is thus a sample analogue of this population identity: it sets the observed average score to zero. The variance of the score is the Fisher information I(θ) = E[s²] = −E[∂s/∂θ], and it controls everything downstream.

The asymptotic theory follows from a Taylor expansion of the score about θ₀. To leading order θ̂ − θ₀ ≈ (1/n ∑ s_i) / I(θ₀), and the numerator is a mean of n zero-mean terms. The central limit theorem then gives √n(θ̂ − θ₀) → N(0, I(θ₀)⁻¹). Equivalently var(θ̂) ≈ 1/[n I(θ₀)], which is precisely the Cramér–Rao lower bound: the MLE is asymptotically efficient, achieving the smallest variance any consistent estimator can. This is why maximum likelihood, not some competitor, is the workhorse of parameter estimation.

Two further virtues matter in practice. First, invariance: if θ̂ maximises the likelihood, then g(θ̂) is the MLE of g(θ) for any function g — one can estimate a rate and read off the mean lifetime as its reciprocal without redoing the fit. Second, maximising the log-likelihood is asymptotically the same as minimising the Kullback–Leibler divergence between the true distribution and the fitted model, so ML is choosing the model closest to reality in an information-theoretic sense.

Common misconceptions. The likelihood is not the probability of the parameter — statements like "θ̂ is 95% likely to be correct" mix up frequentist likelihood with Bayesian posterior probability, which additionally requires a prior. A higher likelihood value does not by itself mean a better model unless the number of parameters is held fixed; comparing models of different dimension requires a penalty (AIC, BIC) or a proper test. And the curvature-based error bars are asymptotic — for small samples or near boundaries they can be wildly optimistic.

Worked examples

Example 1 — Poisson rate from event counts. A detector records the number of decays in five equal one-second windows: 2, 4, 3, 5, 6. Model each count as Poisson(λ) and estimate λ.

f(x; λ) = e^−λ λ^x / x!

Poisson mass function for a count x = 0, 1, 2, … A

ℓ(λ) = ∑_i (−λ + x_i ln λ − ln x_i!) = −nλ + (∑_i x_i) ln λ − ∑_i ln x_i!

Sum the individual log-masses. A

dℓ/dλ = −n + (∑_i x_i)/λ = 0 ⇒ λ̂ = (1/n) ∑_i x_i = x̄

Score equation; the last term is λ-free and drops out. B

d²ℓ/dλ² = −(∑_i x_i)/λ² < 0

Negative for any positive data, confirming a maximum. B

∑ x_i = 2+4+3+5+6 = 20, n = 5 ⇒ λ̂ = 20/5 = 4.0

Insert numbers only now. A

λ̂ = 4.0 events per second

Reading. The best-fit Poisson rate is just the average count. The standard error is √(λ̂/n) = √(4/5) ≈ 0.89 s⁻¹, from I(λ) = n/λ.

Example 2 — Gaussian mean and variance from measurements. Four independent length measurements (metres): 1.70, 1.75, 1.80, 1.65. Model as N(μ, σ²) and estimate both parameters.

ln f(x; μ, σ²) = −½ ln(2πσ²) − (x − μ)² / (2σ²)

Gaussian log-density. A

∂ℓ/∂μ = ∑_i (x_i − μ)/σ² = 0 ⇒ μ̂ = x̄

Score for μ; the common 1/σ² cancels. B

∂ℓ/∂σ² = −n/(2σ²) + ∑_i(x_i − μ)²/(2σ⁴) = 0 ⇒ σ̂² = (1/n) ∑_i(x_i − μ̂)²

Score for σ²; evaluate at μ̂. B

μ̂ = (1.70+1.75+1.80+1.65)/4 = 6.90/4 = 1.725 m

Sample mean. A

deviations: −0.025, +0.025, +0.075, −0.075 m; squares sum to 0.0125 m²

(0.025²+0.025²+0.075²+0.075²) = 0.000625·2 + 0.005625·2 = 0.0125. A

σ̂² = 0.0125/4 = 0.003125 m², σ̂ = 0.0559 m

Divide by n = 4 (the MLE, not the unbiased n−1). A

μ̂ = 1.725 m, σ̂ = 0.056 m

Reading. The mean estimate is the sample average; the MLE variance divides by n and is biased low — the unbiased value would be 0.0125/3 = 0.00417 m² (σ = 0.065 m), the familiar n−1 correction.

Problems

A coin is flipped 10 times, giving 7 heads. Treating each flip as Bernoulli(p), find the MLE of p from the score equation.
Solution
For k successes in n trials, ℓ(p) = k ln p + (n−k) ln(1−p). Then dℓ/dp = k/p − (n−k)/(1−p) = 0. Cross-multiplying: k(1−p) = (n−k)p ⇒ k = np ⇒ p̂ = k/n. With k = 7, n = 10: p̂ = 0.70. Second derivative −k/p² − (n−k)/(1−p)² < 0, a maximum.
Waiting times between five events (hours): 2.1, 1.8, 3.2, 2.5, 2.9. Model as exponential with density f(x; λ) = λe^−λx. Find the MLE of the rate λ.
Solution
ℓ(λ) = ∑(ln λ − λx_i) = n ln λ − λ ∑ x_i. Score: dℓ/dλ = n/λ − ∑ x_i = 0 ⇒ λ̂ = n/∑ x_i = 1/x̄. Sum = 2.1+1.8+3.2+2.5+2.9 = 12.5 h, x̄ = 2.5 h. λ̂ = 5/12.5 = 0.40 per hour (mean waiting time 1/λ̂ = 2.5 h, consistent with invariance).
Prove that the score has zero expectation at the true parameter, E[∂ ln f/∂θ] = 0, stating where regularity is used.
Solution
Normalisation holds for every θ: ∫f(x; θ) dx = 1. Differentiate both sides in θ and interchange ∂/∂θ with the integral (this is the regularity / dominated-convergence step, valid only when the support does not depend on θ): ∫∂f/∂θ dx = 0. Now ∂f/∂θ = f · ∂ ln f/∂θ, so ∫f (∂ ln f/∂θ) dx = E[∂ ln f/∂θ] = 0. This underlies the consistency of the MLE, which matches the sample-average score to this population value.
Three measurements with a known instrument noise σ = 2.0 (units): 5, 7, 9. Assuming N(μ, σ²) with σ fixed, find the MLE of μ and its standard error.
Solution
With σ known, ℓ(μ) = const − ∑(x_i−μ)²/(2σ²). Score ∑(x_i−μ)/σ² = 0 ⇒ μ̂ = x̄ = (5+7+9)/3 = 7.0. Fisher information I = n/σ² = 3/4, so var(μ̂) = σ²/n = 4/3 and the standard error is √(4/3) ≈ 1.15. Result: μ̂ = 7.0 ± 1.2.
(Non-regular.) Data 1.2, 3.7, 4.8, 2.9 are modelled as uniform(0, θ). Show the score equation gives no interior solution and find θ̂ directly; comment on its bias.
Solution
The density is f = 1/θ for 0 ≤ x ≤ θ, else 0. The likelihood is L(θ) = θ⁻ⁿ provided θ ≥ max x_i (otherwise some point falls outside the support and L = 0). Where nonzero, dℓ/dθ = −n/θ, which is never zero — the score equation has no solution. Since L decreases in θ, it is maximised at the smallest permitted value: θ̂ = max x_i = 4.8. This is the boundary MLE, illustrating the support-dependence failure mode. It is biased downward: E[θ̂] = nθ/(n+1) < θ, so an unbiased version multiplies by (n+1)/n = 5/4, giving 6.0.