30. Fisher information and the Cramér–Rao bound

Sufficiency told us which statistics carry the information in a sample. This lesson quantifies how much information there is — and derives a hard limit on how good any unbiased estimator can possibly be.

The score

Start with the derivative of the log-likelihood, called the score:

s(\theta) = \frac{\partial}{\partial\theta}\log f(X \mid \theta)

Note X is random here, so the score is a random variable.

Under mild regularity conditions its mean is zero:

E[s(\theta)] = 0

Why: \int f(x\mid\theta)dx = 1 for every \theta, so differentiating both sides gives \int \frac{\partial f}{\partial\theta}dx = 0. Rewriting \frac{\partial f}{\partial\theta} = f \cdot \frac{\partial \log f}{\partial\theta} turns that integral into exactly E[s(\theta)].

Setting the score to zero is precisely the MLE equation from two lessons ago — the MLE looks for where the observed score vanishes.

Fisher information

Since the score has mean zero, its variance measures how much it fluctuates:

I(\theta) = \operatorname{Var}\big(s(\theta)\big) = E\left[\left(\frac{\partial \log f}{\partial\theta}\right)^2\right]

Under the same regularity conditions there's an equivalent and usually easier form:

I(\theta) = -E\left[\frac{\partial^2}{\partial\theta^2}\log f(X\mid\theta)\right]

The second form is the one to hold in your head. The second derivative measures the curvature of the log-likelihood at its peak.

• Sharply peaked log-likelihood ⟹ large curvature ⟹ high information. Small changes in \theta make the data much less likely, so \theta is pinned down.
• Flat log-likelihood ⟹ low curvature ⟹ little information. Many parameter values explain the data about equally well.

Information is literally the sharpness of the likelihood peak.

It adds up

For n independent observations:

I_n(\theta) = n\,I_1(\theta)

Log-likelihoods add across independent observations, so their curvatures add. Twice the data, twice the information — and this linear growth is where the 1/\sqrt n from §3 ultimately comes from.

The Cramér–Rao lower bound

Here is the payoff. For any unbiased estimator \hat\theta:

\boxed{\;\operatorname{Var}(\hat\theta) \;\ge\; \frac{1}{I_n(\theta)} = \frac{1}{n\,I_1(\theta)}\;}

No unbiased estimator can do better than this. It's not a statement about clever techniques; it's a limit imposed by the data-generating process itself.

An estimator attaining the bound is called efficient. And the general form for a biased estimator of g(\theta):

\operatorname{Var}(\hat\theta) \ge \frac{\big[g'(\theta)\big]^2}{I_n(\theta)}

This is why the MLE is so highly regarded: its asymptotic variance is exactly 1/I_n(\theta), so it attains the bound in the limit — asymptotically efficient, as claimed two lessons ago.

The bound applies to unbiased estimators only. Biased estimators can and do have smaller variance — sometimes smaller MSE than the Cramér–Rao bound itself. That is not a contradiction; it's the bias–variance trade-off from lesson 1 of this module, and it's what ridge regression exploits in §6.

Worked example 1: Bernoulli

\log f(x \mid p) = x\log p + (1-x)\log(1-p)

\frac{\partial \log f}{\partial p} = \frac{x}{p} - \frac{1-x}{1-p}

\frac{\partial^2 \log f}{\partial p^2} = -\frac{x}{p^2} - \frac{1-x}{(1-p)^2}

Take -E[\cdot], using E[X] = p:

I_1(p) = \frac{p}{p^2} + \frac{1-p}{(1-p)^2} = \frac{1}{p} + \frac{1}{1-p} = \frac{1}{p(1-p)}

So the bound for n observations is

\operatorname{Var}(\hat p) \ge \frac{p(1-p)}{n}

And the sample proportion has variance exactly p(1-p)/n. \hat p = \bar X is efficient — no unbiased estimator of p can beat it, at any sample size.

Notice I_1(p) = \frac{1}{p(1-p)} is minimised at p = 0.5 and blows up near 0 or 1. A coin near 50/50 is the hardest case to pin down; a very biased coin reveals itself quickly.

Worked example 2: Normal mean

\log f(x\mid\mu) = -\tfrac12\log(2\pi\sigma^2) - \frac{(x-\mu)^2}{2\sigma^2}

\frac{\partial \log f}{\partial\mu} = \frac{x - \mu}{\sigma^2}, \qquad \frac{\partial^2 \log f}{\partial\mu^2} = -\frac{1}{\sigma^2}

The second derivative is a constant, so taking -E[\cdot] is trivial:

I_1(\mu) = \frac{1}{\sigma^2}, \qquad I_n(\mu) = \frac{n}{\sigma^2}

Bound:

\operatorname{Var}(\hat\mu) \ge \frac{\sigma^2}{n}

And \operatorname{Var}(\bar X) = \sigma^2/n exactly. The sample mean is efficient, which finally justifies its universal use — it isn't just convenient, it is provably optimal among unbiased estimators.

Note also that information is inversely proportional to \sigma^2: noisier data carries less information per observation, exactly as intuition demands.

Where the bound doesn't apply

The regularity conditions matter, and the standard counterexample is the Uniform(0,\theta) again.

The support depends on \theta, so you cannot differentiate under the integral sign — the interchange used to prove E[s] = 0 fails. And indeed the bias-corrected maximum \frac{n+1}{n}X_{(n)} has variance of order 1/n^2, while any Cramér–Rao-style bound would suggest order 1/n.

The estimator beats the "bound" because the bound was never valid here. Any time the support depends on the parameter, the whole framework needs care.

Doing it in Python

Estimate Fisher information numerically as the curvature of the log-likelihood, and compare with theory:

import numpy as np

p_true, n = 0.3, 1
theoretical = 1 / (p_true * (1 - p_true))

# I(p) = -E[d^2/dp^2 log f], computed by averaging over the two outcomes
def d2_logf(x, p, h=1e-5):
    f = lambda q: x * np.log(q) + (1 - x) * np.log(1 - q)
    return (f(p + h) - 2 * f(p) + f(p - h)) / h**2

numerical = -(p_true * d2_logf(1, p_true) + (1 - p_true) * d2_logf(0, p_true))

print("Fisher information for one Bernoulli observation")
print("  theory   1/(p(1-p)) =", round(theoretical, 6))
print("  numeric             =", round(numerical, 6))

Information depends on where you are in the parameter space:

import numpy as np

print(f"{'p':>6} {'I(p) = 1/(p(1-p))':>20} {'CR bound at n=100':>20}")
for p in (0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99):
    info = 1 / (p * (1 - p))
    print(f"{p:>6} {info:>20.4f} {1/(100*info):>20.6f}")

print("\nInformation is LOWEST at p = 0.5 -- a fair coin is the hardest to pin down.")

Confirming the sample proportion attains the bound:

import numpy as np

rng = np.random.default_rng(0)
p, trials = 0.3, 400_000

print(f"{'n':>6} {'Var(p_hat)':>14} {'CR bound':>14} {'ratio':>8}")
for n in (10, 50, 200, 1000):
    phat = rng.binomial(n, p, size=trials) / n
    bound = p * (1 - p) / n
    print(f"{n:>6} {phat.var():>14.8f} {bound:>14.8f} {phat.var()/bound:>8.4f}")

print("\nRatio is 1.00 at every n -- the sample proportion is exactly efficient.")

And the case where the bound doesn't hold, because regularity fails:

import numpy as np

rng = np.random.default_rng(1)
theta, trials = 5.0, 200_000

print("Uniform(0, theta): support depends on theta, so Cramer-Rao does not apply.\n")
print(f"{'n':>6} {'Var of corrected max':>22} {'~1/n scaling would give':>26}")
for n in (10, 20, 40, 80):
    est = (n + 1) / n * rng.uniform(0, theta, size=(trials, n)).max(axis=1)
    print(f"{n:>6} {est.var():>22.6f} {theta**2/(3*n):>26.6f}")

print("\nDoubling n cuts the variance by ~4 (order 1/n^2), not by 2.")

Your turn

1. X \sim \text{Poisson}(\lambda). Find I_1(\lambda) and the Cramér–Rao bound for n observations.

2. Is \bar X efficient for the Poisson mean?

3. If I_1(\theta) = 4 and n = 100, what's the smallest possible standard error of an unbiased estimator?

Solutions

1. For a Poisson, $\log f(x\mid\lambda) = x\log\lambda - \lambda - \log(x!)$.

\frac{\partial \log f}{\partial\lambda} = \frac{x}{\lambda} - 1, \qquad \frac{\partial^2 \log f}{\partial\lambda^2} = -\frac{x}{\lambda^2}

Take -E[\cdot], using E[X] = \lambda:

I_1(\lambda) = \frac{\lambda}{\lambda^2} = \frac{1}{\lambda}

\operatorname{Var}(\hat\lambda) \ge \frac{1}{n/\lambda} = \frac{\lambda}{n}

2. Yes. For a Poisson, \operatorname{Var}(X_i) = \lambda, so

\operatorname{Var}(\bar X) = \frac{\lambda}{n}

which exactly equals the bound. \bar X is efficient — and it's also the MLE and the method-of-moments estimator, all three agreeing.

Note the information 1/\lambda decreases as \lambda grows: rarer events are individually more informative about the rate, because a Poisson with large \lambda is relatively less variable (\sigma/\mu = 1/\sqrt\lambda).

3.

\operatorname{Var}(\hat\theta) \ge \frac{1}{nI_1(\theta)} = \frac{1}{100 \times 4} = 0.0025

\text{SE} \ge \sqrt{0.0025} = 0.05

No unbiased estimator built from these 100 observations can have a standard error below 0.05. To halve it to 0.025 you'd need n = 400 — the familiar \sqrt n tax, here derived from information rather than from the CLT.

Check yourself in code

Verify that the sample proportion attains the Cramér–Rao bound by comparing its simulated variance against p(1-p)/n.

Print exactly this:

n=10 ratio 1.002
n=100 ratio 1.003
n=1000 ratio 1.003
efficient: True

Use default_rng(0), p = 0.3, and 400000 trials per n. Print the ratio of simulated variance to the bound to 3 decimal places. Report efficient: True if every ratio is within 0.01 of 1.

Three decimals is the honest precision here. A variance estimated from 400,000 draws has a relative error near \sqrt{2/400000} \approx 0.2\%, so the ratio is only pinned down to about \pm0.002 — the third decimal. The fourth digit is simulation noise, and it also moves between NumPy versions, since binomial switches algorithm as n grows. A ratio of 1.003 is not evidence of a 0.3% inefficiency; it is evidence of a bound being attained, measured with finite samples.

import numpy as np

rng = np.random.default_rng(0)
p, trials = 0.3, 400_000
ratios = []

for n in (10, 100, 1000):
    phat = rng.binomial(n, p, size=trials) / n
    bound = p * (1 - p) / n
    r = phat.var() / bound
    ratios.append(r)
    print(f"n={n} ratio {r:.3f}")

# Report whether every ratio is within 0.01 of 1.

import numpy as np

rng = np.random.default_rng(0)
p, trials = 0.3, 400_000
ratios = []

for n in (10, 100, 1000):
    phat = rng.binomial(n, p, size=trials) / n
    bound = p * (1 - p) / n
    r = phat.var() / bound
    ratios.append(r)
    print(f"n={n} ratio {r:.3f}")

print("efficient:", all(abs(r - 1) < 0.01 for r in ratios))

---

Fisher information is the curvature of the log-likelihood — how sharply the data pins down the parameter. It adds across independent observations, and its reciprocal is the Cramér–Rao lower bound on the variance of any unbiased estimator. The sample mean and sample proportion attain it exactly; the MLE attains it asymptotically. And when the support depends on the parameter, the bound doesn't apply at all.

Next: turning a point estimate plus its standard error into an interval that states its own uncertainty.