The Neyman–Pearson lemma and likelihood ratio tests · Probability & Statistics

We've been building tests by intuition — reject when \bar X is far from \mu_0. But is that the best test? Among all tests with \alpha = 0.05, which has the most power?

The Neyman–Pearson lemma answers this completely for simple hypotheses, and the answer generalises into the standard recipe for building tests.

Simple hypotheses

A simple hypothesis specifies the distribution completely:

$H_0: \theta = \theta_0 \quad \text{vs} \quad H_1: \theta = \theta_1$

Both are single points — no free parameters left. (A composite hypothesis like \theta > \theta_0 covers a range; that's harder, and comes later in this lesson.)

The lemma

Among all tests with significance level \le \alpha, the most powerful test of H_0 vs H_1 rejects when the likelihood ratio is small:

$\Lambda(x) = \frac{L(\theta_0 \mid x)}{L(\theta_1 \mid x)} \le k$

with k chosen so that P(\Lambda \le k \mid H_0) = \alpha.

The intuition is exactly right: reject when the data is much better explained by H_1 than by H_0. What the lemma adds is that no other test can do better — this rule is optimal, not merely reasonable.

Why it's optimal

The proof is a neat exchange argument. Think of building a rejection region by choosing which data points to include, subject to a budget: the total probability under H_0 can't exceed \alpha.

Each point x costs f(x \mid \theta_0) from the budget and returns f(x \mid \theta_1) in power. To maximise the return, take the points with the best ratio of return to cost — largest f(x\mid\theta_1)/f(x\mid\theta_0), i.e. smallest \Lambda — until the budget is spent.

That's the greedy solution to a knapsack problem, and here it's exactly optimal because the items are infinitesimally divisible.

Worked example: Normal mean

X_1,\dots,X_n \sim N(\mu, \sigma^2) with \sigma known. Test H_0: \mu = 0 vs H_1: \mu = 1.

$\Lambda = \frac{\exp\left(-\frac{1}{2\sigma^2}\sum x_i^2\right)}{\exp\left(-\frac{1}{2\sigma^2}\sum(x_i - 1)^2\right)}$

Take logs:

$\log\Lambda = -\frac{1}{2\sigma^2}\left[\sum x_i^2 - \sum(x_i-1)^2\right] = -\frac{1}{2\sigma^2}\left[2\sum x_i - n\right]$

So \log \Lambda is a decreasing linear function of \sum x_i. Therefore

$\Lambda \le k \iff \sum x_i \ge c \iff \bar X \ge c'$

The optimal test rejects for large \bar X — exactly the test we'd been using by intuition. Now we know it's not just reasonable, it's the most powerful test that exists at that \alpha.

Choose c' from the null distribution: c' = z_{1-\alpha}\sigma/\sqrt n.

Notice the cutoff doesn't depend on \mu_1 = 1 at all. Any \mu_1 > 0 gives the same test — which is what makes the next step possible.

Uniformly most powerful tests

A test that is most powerful for every alternative in a composite H_1 is uniformly most powerful (UMP).

The example above is UMP for H_0: \mu = 0 vs H_1: \mu > 0, because the same rejection region is optimal against every \mu_1 > 0 simultaneously.

This is why one-sided tests are well behaved and two-sided ones are not. For H_1: \mu \ne 0, the optimal test against \mu_1 = +1 rejects for large \bar X, while against \mu_1 = -1 it rejects for small \bar X. No single region is optimal for both, so no UMP test exists for a two-sided alternative. The symmetric two-tailed test is a sensible compromise, not an optimum.

UMP tests exist whenever the family has a monotone likelihood ratio — the ratio L(\theta_1)/L(\theta_0) is monotone in some statistic T. That covers most exponential families.

The generalised likelihood ratio test

Real problems have nuisance parameters and composite hypotheses on both sides. The general recipe: maximise the likelihood under each hypothesis and compare.

$\Lambda = \frac{\sup_{\theta \in \Theta_0} L(\theta)}{\sup_{\theta \in \Theta} L(\theta)}$

The numerator is the best the null can do; the denominator is the best anything can do. So 0 \le \Lambda \le 1, and small \Lambda means the restriction to H_0 costs a lot — evidence against it.

Wilks' theorem

The distribution of \Lambda is usually intractable, but its asymptotic behaviour is not. Under H_0, with mild regularity conditions:

$-2\log\Lambda \;\xrightarrow{d}\; \chi^2_r$

where r is the difference in the number of free parameters between the full model and the null.

That single result is enormous in practice. It means you can:

Fit the model under H_0, record the maximised log-likelihood \ell_0.
Fit the full model, record \ell_1.
Compute -2(\ell_0 - \ell_1) and compare to \chi^2_r.

This works for essentially any model you can fit by maximum likelihood — logistic regression, GLMs, mixed models. It's the machinery behind the deviance tests in §6 and the model comparisons throughout modern statistics.

Worked example: the LRT for a Normal variance

H_0: \sigma^2 = \sigma_0^2 vs H_1: \sigma^2 \ne \sigma_0^2, with \mu unknown.

Under H_1 the MLEs are \hat\mu = \bar x and $\hat\sigma^2 = \frac1n\sum(x_i - \bar x)^2$. Under H_0, \hat\mu = \bar x still, but \sigma^2 is fixed.

Writing W = \frac{n\hat\sigma^2}{\sigma_0^2}, the ratio works out to

$-2\log\Lambda = W - n - n\log\frac{W}{n}$

One free parameter is constrained (\sigma^2), so r = 1 and

$-2\log\Lambda \;\xrightarrow{d}\; \chi^2_1$

The exact test here uses W \sim \chi^2_{n-1} directly, so we can check how good the asymptotic approximation is — which the code below does.

Doing it in Python

Confirming the Neyman–Pearson test really is most powerful, by comparing it against a plausible alternative test at the same \alpha:

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(0)
n, sigma, mu1, alpha, trials = 20, 1.0, 0.6, 0.05, 200_000

# Neyman-Pearson: reject for large xbar
crit_np = norm.ppf(1 - alpha) * sigma / np.sqrt(n)

# A reasonable competitor: reject for a large MEDIAN
null = rng.normal(0, sigma, size=(trials, n))
crit_med = np.quantile(np.median(null, axis=1), 1 - alpha)

alt = rng.normal(mu1, sigma, size=(trials, n))
power_np = (alt.mean(axis=1) > crit_np).mean()
power_med = (np.median(alt, axis=1) > crit_med).mean()

print(f"both tests calibrated to alpha = {alpha}")
print(f"  Neyman-Pearson (mean)  power: {power_np:.4f}")
print(f"  competitor (median)    power: {power_med:.4f}")
print("\nNo test at this alpha can beat the first -- that's the lemma.")

Wilks' theorem, checked against its asymptotic claim:

import numpy as np
from scipy.stats import chi2, kstest

rng = np.random.default_rng(1)
n, sigma0, trials = 50, 2.0, 100_000

x = rng.normal(5.0, sigma0, size=(trials, n))          # H0 is TRUE
sigma_hat2 = x.var(axis=1, ddof=0)
W = n * sigma_hat2 / sigma0**2

stat = W - n - n * np.log(W / n)

print("Wilks: -2 log Lambda should be chi2 with 1 df")
print("  mean:", round(stat.mean(), 4), " (chi2_1 mean is 1)")
print("  var :", round(stat.var(), 4), " (chi2_1 var is 2)")
for q in (0.5, 0.9, 0.95, 0.99):
    print(f"  q={q}: simulated {np.quantile(stat, q):.4f}   chi2_1 {chi2.ppf(q, 1):.4f}")

Building a likelihood ratio test from scratch for a real comparison — two Poisson rates:

import numpy as np
from scipy.stats import chi2

rng = np.random.default_rng(2)

def lrt_two_poisson(x, y):
    """H0: both samples share one rate. H1: separate rates."""
    n, m = len(x), len(y)
    lam_pooled = (x.sum() + y.sum()) / (n + m)
    lam_x, lam_y = x.mean(), y.mean()

    def ll(data, lam):
        return -len(data) * lam + data.sum() * np.log(lam) if lam > 0 else -np.inf

    ll0 = ll(x, lam_pooled) + ll(y, lam_pooled)
    ll1 = ll(x, lam_x) + ll(y, lam_y)
    stat = -2 * (ll0 - ll1)
    return stat, 1 - chi2.cdf(stat, df=1)      # 2 free params vs 1, so r = 1

# Same rate: should not reject
a = rng.poisson(3.0, 200)
b = rng.poisson(3.0, 200)
s, p = lrt_two_poisson(a, b)
print(f"same rate     : stat {s:8.4f}  p {p:.4f}")

# Different rates: should reject
c = rng.poisson(3.0, 200)
d = rng.poisson(4.2, 200)
s, p = lrt_two_poisson(c, d)
print(f"different rates: stat {s:8.4f}  p {p:.6f}")

Your turn

1. Why does no UMP test exist for H_0: \mu = 0 vs H_1: \mu \ne 0?

2. For the LRT, why is \Lambda \le 1 always?

3. You compare a model with 3 parameters against one with 7. What are the degrees of freedom for Wilks' theorem?

Solutions

1. Because the optimal rejection region points in different directions for different alternatives.

Against \mu_1 = +1, Neyman–Pearson says reject for large \bar X. Against \mu_1 = -1, the same argument says reject for small \bar X. These are different regions, and neither is most powerful against the other's alternative.

Since UMP requires one region to be simultaneously optimal against every alternative in H_1, and no such region exists here, there is no UMP test.

The usual two-tailed test splits \alpha evenly and is a good compromise — it's UMP among unbiased tests, a weaker but still meaningful optimality.

2. Because the numerator maximises over a subset of the space the denominator maximises over:

$\Theta_0 \subseteq \Theta \implies \sup_{\theta\in\Theta_0} L(\theta) \le \sup_{\theta\in\Theta} L(\theta)$

A constrained maximum can never exceed an unconstrained one. So the ratio is at most 1, hitting 1 exactly when the unrestricted MLE happens to satisfy H_0 anyway — in which case the null costs nothing and there's no evidence against it.

3. r = 7 - 3 = 4 degrees of freedom.

The statistic -2\log\Lambda is compared against \chi^2_4. Intuitively, each constrained parameter contributes one degree of freedom, because each one is free to move in the full model and pinned in the reduced one.

An important caveat: this requires the models to be nested (the smaller is a special case of the larger) and the null parameters to lie in the interior of the space. Testing whether a variance component is zero puts the null on a boundary, and the \chi^2_r approximation fails there.

Check yourself in code

Verify Wilks' theorem: simulate under a true null and confirm -2\log\Lambda follows a \chi^2_1 distribution.

Print exactly this:

mean 1.0395
var 2.14
q95 simulated 4.0171
q95 chi2_1 3.8415
matches: True

Use default_rng(1), n = 50, \sigma_0 = 2, true mean 5, and 100000 trials. Round every value to 4 decimal places. Report matches: True if the simulated 95th percentile is within 0.2 of the \chi^2_1 value.

Note the gap: at n = 50 the simulated 95th percentile is 4.02 against the asymptotic 3.84. Wilks' theorem is a limit, and at this sample size the approximation is close but not exact — using it here would reject slightly too often, giving a real Type I rate nearer 6% than 5%. Raise n to 500 and the two agree to two decimal places.

import numpy as np
from scipy.stats import chi2

rng = np.random.default_rng(1)
n, sigma0, trials = 50, 2.0, 100_000

x = rng.normal(5.0, sigma0, size=(trials, n))
W = n * x.var(axis=1, ddof=0) / sigma0**2
stat = W - n - n * np.log(W / n)

print("mean", round(stat.mean(), 4))

# Print the variance, the simulated 95th percentile, the chi2_1 95th
# percentile, and whether they agree to within 0.2.

import numpy as np
from scipy.stats import chi2

rng = np.random.default_rng(1)
n, sigma0, trials = 50, 2.0, 100_000

x = rng.normal(5.0, sigma0, size=(trials, n))
W = n * x.var(axis=1, ddof=0) / sigma0**2
stat = W - n - n * np.log(W / n)

print("mean", round(stat.mean(), 4))
print("var", round(stat.var(), 4))

q95 = np.quantile(stat, 0.95)
theory = chi2.ppf(0.95, 1)
print("q95 simulated", round(q95, 4))
print("q95 chi2_1", round(theory, 4))
print("matches:", bool(abs(q95 - theory) < 0.2))

For simple hypotheses, the likelihood ratio test is provably the most powerful test at any given \alpha — that's the Neyman–Pearson lemma. When the same region stays optimal across a whole composite alternative you get a UMP test, which is why one-sided tests behave better than two-sided ones. And the generalised LRT, with Wilks' \chi^2 approximation, extends the idea to essentially any pair of nested models you can fit.

Next: the workhorse application — testing means when \sigma is unknown.

35. The Neyman–Pearson lemma and likelihood ratio tests