Parameter Uncertainties and the Cramer-Rao Bound

Derivation

Parameter Uncertainties and the Cramer-Rao Bound

D-122 Home PU-106 Threads chance Depends on The Method of Maximum Likelihood, Linear Least Squares and the Normal Equations

Statement

For a statistical model with log-likelihood ℓ(θ) depending on a parameter vector θ = (θ₁, …, θ_p), we derive that (i) the maximum-likelihood estimator θ̂ has, asymptotically, covariance equal to the inverse of the Fisher information matrix I(θ), where I_jk = −E[∂²ℓ/∂θ_j∂θ_k] is the expected curvature of the log-likelihood; and (ii) for any unbiased estimator the covariance is bounded below by I⁻¹ in the matrix (Löwner) sense — the Cramér–Rao bound cov(θ̂) ⪰ I(θ)⁻¹.

Why it matters

Every measured "value ± error" quoted from a fit rests on this result. The Fisher matrix converts the shape of the likelihood surface near its peak into the covariance ellipsoid of the fitted parameters, so a physicist can read off parameter uncertainties and their correlations directly from the second derivatives of ℓ — no repeated experiments required.

The Cramér–Rao bound sets the best variance any unbiased method can achieve, so it defines the target against which a real estimator is judged "efficient". It underlies error forecasting for future experiments (Fisher forecasting in cosmology and particle physics) and the design of measurements that minimise a chosen parameter's variance.

Assumptions

The likelihood L(θ) = p(x | θ) is a correctly specified probability density for the data x.If the model is wrong, curvature still gives an error bar but it estimates the spread about the wrong value; the bound loses its meaning and the "sandwich" covariance replaces I⁻¹.

The estimator considered for the bound is unbiased: E[θ̂] = θ for all θ.Biased estimators can beat I⁻¹ (e.g. shrinkage); the bound then applies to the bias-corrected form (∂E[θ̂]/∂θ) I⁻¹ (∂E[θ̂]/∂θ)^T.

Regularity: the support of p(x|θ) does not depend on θ, and differentiation and integration over x may be interchanged.If dropped (e.g. a uniform density on [0, θ]) the score identity E[∂ℓ/∂θ] = 0 fails and the Cramér–Rao bound does not hold — such models can converge at 1/n instead of 1/√n.

Local quadratic (asymptotic normal) behaviour: near the peak ℓ is well approximated by its second-order Taylor expansion.If dropped — small samples, boundaries, degeneracies — the covariance-from-curvature identification is only leading order and the true confidence regions are non-elliptical.

Derivation

s(θ) ≡ ∇_θ ℓ(θ) = ∇_θ log p(x | θ), s_j = ∂ℓ/∂θ_j

Define the score vector, the gradient of the log-likelihood. This is the object whose zero locates the MLE. A

∫ p(x | θ) dx = 1 ⟹ ∂/∂θ_j ∫ p dx = ∫ ∂p/∂θ_j dx = 0

Normalisation holds for every θ; differentiate both sides and interchange derivative and integral (regularity assumption). B

∫ ∂p/∂θ_j dx = ∫ (∂ log p/∂θ_j) p dx = E[s_j] = 0

Use the log-derivative identity ∂p/∂θ = p · ∂log p/∂θ. The score has zero mean — the pivotal fact of the whole derivation. B

∂/∂θ_k ∫ s_j p dx = ∫ (∂s_j/∂θ_k) p dx + ∫ s_j (∂log p/∂θ_k) p dx = 0

Differentiate the identity E[s_j] = 0 a second time (again interchanging), and rewrite ∂p/∂θ_k = p s_k. C

E[s_j s_k] = −E[∂s_j/∂θ_k] = −E[∂²ℓ/∂θ_j∂θ_k] ≡ I_jk

Rearranging step 4 gives the two equivalent faces of the Fisher information: variance of the score equals expected negative curvature (the Hessian) of ℓ. Since E[s]=0, I = cov(s) is symmetric positive semi-definite. C

0 = s(θ̂) ≈ s(θ) + H(θ)(θ̂ − θ), H_jk = ∂²ℓ/∂θ_j∂θ_k

The MLE satisfies s(θ̂)=0 (prior result: maximum-likelihood estimation). Taylor-expand the score to first order about the true θ; legal in the local-quadratic regime. B

θ̂ − θ ≈ −H(θ)⁻¹ s(θ)

Solve the linearised score equation for the estimator error. This is the exact analogue of solving the normal equations in linear least squares (prior result), with −H playing the role of the design matrix A^TA. B

cov(θ̂) ≈ H⁻¹ cov(s) H^−T = H⁻¹ I H⁻¹

Propagate covariance through the linear map θ̂−θ = −H⁻¹s using cov(Ms)=M cov(s) M^T and cov(s)=I from step 5. B

E[H] = −I ⟹ cov(θ̂) ≈ I⁻¹ I I⁻¹ = I(θ)⁻¹

In the asymptotic regime the observed curvature concentrates on its expectation −I (law of large numbers), so the sandwich collapses to the inverse Fisher matrix. This is the asymptotic MLE covariance. B

cov [ θ̂ ⊕ s ] = [ Σ C ; C^T I ] ⪰ 0, where C_jk = cov(θ̂_j, s_k)

Now derive the bound for a general unbiased estimator, not just the MLE. Any covariance matrix of the stacked vector (θ̂, s) is positive semi-definite. Write Σ = cov(θ̂). C

C_jk = E[θ̂_j s_k] = ∂/∂θ_k E[θ̂_j] = ∂θ_j/∂θ_k = δ_jk

Unbiasedness gives E[θ̂_j] = θ_j; differentiate under the integral, using ∂p/∂θ_k=p s_k and E[s]=0. The cross-covariance is the identity. C

[ Σ 𝟙 ; 𝟙 I ] ⪰ 0 ⟹ Σ − 𝟙 · I⁻¹ · 𝟙 = Σ − I⁻¹ ⪰ 0

Apply the Schur-complement condition for positive semi-definiteness of a block matrix (with I ≻ 0). The complement of the I-block must be PSD. C

cov(θ̂) ⪰ I(θ)⁻¹

Rearrange: the covariance of any unbiased estimator exceeds the inverse Fisher information in the Löwner order. The MLE (step 9) attains it asymptotically, hence is asymptotically efficient. C

Result

I_jk(θ) = −E[ ∂²ℓ/∂θ_j∂θ_k ] = E[ s_j s_k ] ⟹ cov(θ̂) ⪰ I(θ)⁻¹, cov(θ̂_MLE) → I(θ)⁻¹

Reading. The Fisher matrix is the average sharpness of the likelihood peak — large curvature (a narrow, steep peak) means the data pin the parameters down tightly, so the inverse (the covariance) is small. Off-diagonal elements encode parameter degeneracy: a shallow direction of the likelihood surface (near-singular I) becomes a long axis of the error ellipse. No unbiased estimator can do better than I⁻¹, and maximum likelihood reaches this floor for large samples. The diagonal σ_j = √(I⁻¹)_jj is the marginal 1σ error on θ_j (with the others floating); 1/√I_jj is the tighter error with the others held fixed.

Units check. ℓ is dimensionless (log of a probability), so ∂²ℓ/∂θ_j∂θ_k carries units of [θ_jθ_k]⁻¹. Thus I_jk has units [θ_jθ_k]⁻¹ and its inverse (I⁻¹)_jk has units [θ_jθ_k] — exactly a covariance, with √(I⁻¹)_jj in units of θ_j as an error bar must be.

Limiting cases

Single parameter: I = −E[ℓ″] is a scalar and var(θ̂) ≥ 1/I, σ_θ ≈ 1/√I.
Independent measurements: Fisher information adds, I_total = Σ_i I_i = n I₁, so σ ∝ 1/√n — the canonical statistical scaling.
Gaussian data with known covariance C, mean μ(θ): I_jk = (∂μ/∂θ_j)^T C⁻¹ (∂μ/∂θ_k) — the least-squares curvature, recovering the normal-equations covariance (A^TC⁻¹A)⁻¹.
Orthogonal parameters: I diagonal ⟹ errors uncorrelated, marginal and conditional errors coincide.
Estimating a mean, N(θ, σ²): I = n/σ², bound var ≥ σ²/n, attained exactly by the sample mean for all n.

Breaks when

Boundary / edge-supported models. If the support depends on θ (e.g. Uniform[0,θ]), the regularity interchange of step 2 fails, E[s]≠0, and the Cramér–Rao bound is simply false — the max-order statistic estimates θ with variance ∝ 1/n², far below any I⁻¹ ∝ 1/n.
Singular or near-singular Fisher matrix (degeneracy). When two parameters affect the data almost identically, I is nearly non-invertible; I⁻¹ blows up and the local quadratic picture is invalid. The bound is formally infinite along the flat direction and the reported covariance is dominated by numerical conditioning.
Small samples / non-quadratic peak. The Taylor expansion of step 6 and the concentration H→−I of step 9 require large n; for finite samples the true confidence region is non-elliptical and asymmetric, and the MLE is biased at order 1/n, so I⁻¹ understates the true spread.
Misspecified model. If p(x|θ) is not the true data law, the Hessian and score covariance differ, E[H] ≠ −cov(s), and the correct covariance is the robust sandwich H⁻¹ cov(s) H⁻¹, not I⁻¹.

Failure modes

Inverting element-wise: taking σ_j = 1/√I_jj instead of √(I⁻¹)_jj. The first is the error with other parameters fixed; the correct marginal error requires inverting the full matrix and is always ≥ it.
Sign error: forgetting the minus sign, using +E[ℓ″]. Since the Hessian at a maximum is negative-definite, this yields a negative "information" and imaginary error bars.
Confusing observed and expected information: using the Hessian at θ̂ (observed) versus its expectation I (expected). Both are valid asymptotically but must not be mixed within one covariance.
Applying the bound to a biased estimator and concluding it "violates" Cramér–Rao. The bound only constrains unbiased estimators (or the bias-corrected form).
Treating I⁻¹ as exact for small n, quoting symmetric ± errors when the likelihood is skewed — profiling the likelihood is required instead.
Forgetting information adds: scaling errors as 1/n instead of 1/√n by inverting summed variances rather than summed informations.

Discussion

The deep content of this result is that a single, local, differential quantity — the curvature of ℓ at its peak — controls the global statistical precision of a measurement. The two faces of I, "variance of the score" and "expected negative Hessian", are forced to agree by nothing more than the normalisation of probability differentiated twice (steps 3–5). This is why one can forecast an experiment's error bars from the model alone, before any data are taken: the Fisher matrix depends only on p(x|θ).

Geometrically, I(θ) is a metric on the manifold of probability distributions — the Fisher–Rao metric of information geometry — measuring how distinguishable nearby parameter values are. Directions in which the data barely change the distribution have small information and large uncertainty; the error ellipsoid is literally the unit ball of the inverse metric. This connects parameter estimation to the natural gradient in optimisation and to the KL divergence, whose local expansion is D_KL(θ ‖ θ+dθ) ≈ ½ dθ^T I dθ.

The tie to the assumed prior results is exact. Maximum likelihood provides the estimating equation s(θ̂)=0; linearising it (steps 6–7) is structurally identical to solving the least-squares normal equations, with −H replacing A^TA. For Gaussian data the two coincide precisely: I = A^TC⁻¹A and the ML covariance is the familiar (A^TC⁻¹A)⁻¹. Cramér–Rao then tells us this least-squares covariance is not just convenient but optimal among unbiased estimators when the noise is Gaussian.

The bound is a strict inequality in the Löwner (matrix) sense, which is stronger than bounding each variance separately: it constrains every linear combination a^Tθ̂ simultaneously via var(a^Tθ̂) ≥ a^TI⁻¹a. Equality throughout requires the score to be an affine function of the estimator, s(θ) = I(θ)(θ̂ − θ), which holds identically only for exponential-family models in their natural parametrisation — the precise class for which an efficient unbiased estimator exists at finite n. Outside that class efficiency is only asymptotic, and higher-order (Bhattacharyya) bounds tighten the CR floor. Common misconceptions: that I⁻¹ is an exact covariance for any sample size (it is a large-n limit); that a smaller variance than I⁻¹ is impossible (biased estimators routinely achieve it); and that off-diagonal I_jk=0 means the parameters are physically independent (it means the information factorises locally, nothing more).

Worked examples

Example 1 — Exponential lifetimes. n decays with i.i.d. times t_i ~ (1/τ) e^−t/τ; find the error on the mean lifetime τ.

Model: rate parameter unknown, estimate τ from a sample. A

ℓ(τ) = Σ_i [ −log τ − t_i/τ ] = −n log τ − (Σt_i)/τ

Sum the log-density over independent events. A

∂ℓ/∂τ = −n/τ + (Σt_i)/τ², ∂²ℓ/∂τ² = n/τ² − 2(Σt_i)/τ³

First derivative set to zero gives τ̂ = (Σt_i)/n = t̄; second derivative for curvature. B

I(τ) = −E[∂²ℓ/∂τ²] = −n/τ² + 2 n E[t]/τ³ = −n/τ² + 2n/τ² = n/τ²

Use E[t_i] = τ. Information grows linearly with sample size. B

σ_τ = √(I⁻¹) = τ/√n. Numbers: τ̂ = 2.30 μs, n = 400 ⟹ σ_τ = 2.30/√400 = 0.115 μs

Invert scalar information and insert the fitted value for τ. A

τ = 2.30 ± 0.12 μs (σ_τ = τ/√n)

Reading. The exponential model is an exponential family in 1/τ, so the sample mean is efficient and the bound is attained even at finite n; the 5% fractional error follows purely from 1/√400.

Example 2 — Straight-line fit, slope–intercept correlation. y_i = a + b x_i + ε_i, ε_i ~ N(0, σ²), known σ; find cov(â, b̂).

Gaussian noise ⟹ Fisher matrix is the least-squares curvature. A

ℓ = −(1/2σ²) Σ_i (y_i − a − b x_i)² + const

Log of the Gaussian likelihood; drop θ-independent normalisation. A

I = (1/σ²) [ n Σx_i ; Σx_i Σx_i² ]

Second derivatives: −∂²ℓ/∂a² = n/σ², −∂²ℓ/∂a∂b = Σx_i/σ², −∂²ℓ/∂b² = Σx_i²/σ². B

cov = I⁻¹ = (σ²/Δ) [ Σx_i² −Σx_i ; −Σx_i n ], Δ = n Σx_i² − (Σx_i)²

Invert the 2×2 matrix; Δ = n²·var(x) up to sign. B

n = 10, σ = 0.20, x = 0,1,…,9 ⟹ Σx = 45, Σx² = 285, Δ = 10·285 − 45² = 825

Plug in the design points. A

σ_a² = 0.04·285/825 = 0.01382, σ_b² = 0.04·10/825 = 4.85×10⁻⁴, cov(a,b) = −0.04·45/825 = −2.18×10⁻³

Read off matrix entries times σ²=0.04. A

σ_a = 0.118, σ_b = 0.0220, ρ_ab = cov/(σ_aσ_b) = −0.84

Reading. Intercept and slope are strongly anti-correlated because the x-values are all positive — raising the slope must lower the fitted intercept to keep the line through the data. Centring x (so Σx_i=0) diagonalises I and removes the correlation, the standard trick for decoupling fit parameters.

Problems

(A) Poisson counts. A detector records N ~ Poisson(μ) counts. Show I(μ) = 1/μ for a single measurement and hence give σ_μ for one exposure yielding N = 400.
Solution
ℓ = N log μ − μ − log N!. Then ∂ℓ/∂μ = N/μ − 1, ∂²ℓ/∂μ² = −N/μ². So I(μ) = −E[−N/μ²] = E[N]/μ² = μ/μ² = 1/μ. Thus σ_μ = √μ ≈ √N = √400 = 20, i.e. μ = 400 ± 20 — the familiar √N counting error, here derived as the Cramér–Rao floor.
(B) Gaussian variance. For n i.i.d. samples from N(0, σ²) with unknown σ²=v, find I(v) and the CR bound on var(v̂).
Solution
ℓ = −(n/2)log v − (1/2v)Σx_i². ∂ℓ/∂v = −n/2v + (Σx²)/2v²; ∂²ℓ/∂v² = n/2v² − (Σx²)/v³. Taking E[Σx²]=nv: I(v) = −(n/2v² − nv/v³) = −(n/2v² − n/v²) = n/2v². Bound: var(v̂) ≥ 2v²/n. The MLE v̂ = (1/n)Σx² has variance exactly 2v²/n, attaining the bound; the unbiased 1/(n−1) estimator has slightly larger variance, illustrating the bias–variance interplay.
(C) Angular resolution. A photon direction is measured as φ_i ~ N(φ, s²), n independent hits, known per-hit resolution s = 3.0°, n = 36. Give the error on the reconstructed direction and state which assumption guarantees the bound is attained.
Solution
Single-parameter Gaussian mean: I(φ) = n/s², so σ_φ = s/√n = 3.0/√36 = 0.50°. The mean is an exponential-family sufficient statistic, so the bound is attained exactly (not just asymptotically) — the sample mean is the efficient estimator for all n.
(B) Marginal vs conditional error. A two-parameter fit yields I = [[4, 3],[3, 4]] (units of parameter⁻²). Compute the marginal error on θ₁ and the conditional error (with θ₂ fixed), and comment.
Solution
det I = 16 − 9 = 7. I⁻¹ = (1/7)[[4, −3],[−3, 4]], so (I⁻¹)₁₁ = 4/7 = 0.571 ⟹ marginal σ₁ = 0.756. Conditional error (other parameter fixed) is 1/√I₁₁ = 1/√4 = 0.500. The marginal error is larger (0.756 > 0.500) because the correlation ρ = −3/4 inflates the uncertainty when θ₂ is allowed to float — the general rule that marginalising over nuisance parameters never decreases an error.
(C) Reparametrisation. Show the Fisher information transforms as I(η) = (dθ/dη)² I(θ) under a smooth 1-D reparametrisation θ = θ(η), and verify that the CR bound is preserved. Apply to θ = τ (lifetime) and η = λ = 1/τ (rate) for the exponential model of Worked Example 1.
Solution
By the chain rule s_η = ∂ℓ/∂η = (dθ/dη)(∂ℓ/∂θ) = (dθ/dη)s_θ. Since I = var(s), I(η) = (dθ/dη)² I(θ). The bound: var(η̂) ≥ 1/I(η) = (dη/dθ)²/I(θ) = (dη/dθ)² var(θ̂), which is exactly the first-order error-propagation rule — so the bound is covariant, i.e. consistent under change of variables. For the exponential model I(τ) = n/τ²; with λ = 1/τ, dτ/dλ = −1/λ² = −τ², so I(λ) = τ⁴·(n/τ²) = n τ² = n/λ². Hence σ_λ = λ/√n — the same fractional error 1/√n as for τ, as it must be for a reciprocal.