Likelihood Ratio Tests

Suppose you fit two nested models to the same data. One is a simple model with kk free parameters. The other is a richer model with k+rk+r. Suppose the richer model fits noticeably better. Is the improvement real? Or could it plausibly be explained by chance alone under the simpler model? This is the single most common question in applied statistics. It arises in genetics (does adding an interaction term help?). It arises in econometrics (does a restricted regression suffice?). It also arises in machine learning (is a mixture-of-two-Gaussians better than one Gaussian?). The likelihood ratio test (LRT) answers it with a universal recipe. Compare the maximum likelihood attainable under the small model to the maximum likelihood attainable under the big model. Then calibrate the resulting ratio against a chi-squared distribution. The idea traces back to the Neyman-Pearson Lemma. It shows that for testing one simple hypothesis against another, the likelihood ratio is the optimal test statistic. The generalization to composite hypotheses is due to Samuel Wilks in 1938. Composite hypotheses pin down a parameter only up to a lower-dimensional constraint, rather than to a single point. This generalization is one of the great unifying results of statistical theory. It explains why so many "different" classical tests are really the same idea in disguise. These include the tt-test, the FF-test for nested linear models, and the chi-squared test of independence in contingency tables. Each is a special or asymptotic case of the generalized likelihood ratio statistic. The big picture is this: maximum likelihood estimation gives us, for any model, a well-defined "best fit" score. Nested models can be compared by how much that score improves. We check this improvement when we allow the extra parameters to move. Wilks's theorem says that under mild regularity conditions, something remarkable happens. Twice the log of this improvement converges, as the sample size grows, to a chi-squared distribution. Its degrees of freedom equal exactly the number of extra parameters, rr. This is remarkable: the same asymptotic distribution applies essentially universally. This holds regardless of the specific parametric family. It requires only that the null hypothesis lies in the interior of the parameter space. It also requires that the usual smoothness conditions for Maximum Likelihood Estimation hold. This article builds the generalized likelihood ratio statistic from first principles. It states and proves Wilks's theorem via a careful Taylor expansion argument. It works through two genuinely asymmetric numerical examples: an exponential rate test and a multinomial goodness-of-fit test. Finally, it shows exactly what breaks — and why — when the regularity conditions fail.

Key takeaways
  • The likelihood ratio statistic compares the best achievable likelihood under H0H_0 to the best achievable likelihood under the full (unrestricted) model.
  • Wilks's theorem: under regularity conditions, 2log⁡Λ→χr22\log \Lambda \to \chi^2_r in distribution under H0H_0, where rr is the number of constraints imposed by H0H_0.
  • The test rejects H0H_0 for large values of −2log⁡Λ-2\log\Lambda; it reduces to familiar tt, zz, and FF tests in classical exponential-family settings.
  • The proof rests on a second-order Taylor expansion of the log-likelihood around the MLE, using consistency and asymptotic normality of the MLE.
  • Regularity matters: if the true parameter sits on the boundary of the parameter space, or the model is non-identifiable, the chi-squared limit fails outright.
  • LRTs underlie model-selection tools (AIC, deviance tests in GLMs) and generalize the classical tests you already know into one master framework.
  • Key idea: the log-likelihood surface's curvature near its maximum encodes exactly how much evidence the data carry against a constraint, and that curvature is universally chi-squared in the large-sample limit.

Where this lives today

The likelihood ratio test is the backbone of deviance-based model comparison in generalized linear models: the "deviance" reported by any GLM software (logistic regression, Poisson regression, etc.) is precisely G=−2log⁡ΛG = -2\log\Lambda comparing a fitted model to the saturated model, and nested-model comparisons (does adding a covariate help?) are literally likelihood ratio tests, referenced against χr2\chi^2_r. In genetics, LOD scores in linkage analysis are rescaled likelihood ratio statistics, and genome-wide association studies routinely use LRTs to compare additive vs. dominance models. In econometrics, LR, Wald, and Lagrange-multiplier (score) tests form a standard trio taught alongside the classical FF-test for linear restrictions, which is the exact (non-asymptotic) analogue of the LRT in Gaussian linear models. In phylogenetics and evolutionary biology, likelihood ratio tests compare nested substitution models of DNA evolution to decide, for instance, whether a molecular clock assumption is justified. In signal processing and radar detection, the Neyman–Pearson-optimal detector for known signals in Gaussian noise is literally a likelihood ratio test, and its generalized version underlies detection theory when signal parameters are unknown (GLRT detectors). In machine learning, likelihood ratio ideas underlie information criteria (AIC is derived from an asymptotic bias-correction to the likelihood ratio) and permeate work on model selection, anomaly detection, and hypothesis testing for mixture and latent-variable models — precisely the boundary cases where the classical theorem, as the counterexample above shows, must be modified. The theorem's asymptotic universality is also the conceptual bridge to Bayesian Inference: An Introduction, where Bayes factors serve an analogous "evidence comparison" role but replace the maximized likelihood with an integrated likelihood over a prior.

Ask the tutor: How does the likelihood ratio test connect to Bayesian model comparison via Bayes factors?

DefinitionLikelihood ratio test statistic

Let X1,…,XnX_1,\dots,X_n be i.i.d. observations with joint density/mass function f(x;θ)f(x;\theta), θ∈Θ⊆Rk\theta \in \Theta \subseteq \mathbb{R}^k. Let L(θ)=∏i=1nf(Xi;θ)L(\theta) = \prod_{i=1}^n f(X_i;\theta) denote the likelihood function and ℓ(θ)=log⁡L(θ)\ell(\theta) = \log L(\theta) the log-likelihood. Consider testing

H0: θ∈Θ0vs.H1: θ∈Θ∖Θ0,H_0:\ \theta \in \Theta_0 \qquad \text{vs.} \qquad H_1:\ \theta \in \Theta \setminus \Theta_0,

where Θ0⊂Θ\Theta_0 \subset \Theta is a restricted (typically lower-dimensional) subset of the parameter space. The generalized likelihood ratio is

Λ=Λ(X1,…,Xn)=sup⁡θ∈Θ0L(θ)sup⁡θ∈ΘL(θ).\Lambda = \Lambda(X_1,\dots,X_n) = \frac{\displaystyle\sup_{\theta \in \Theta_0} L(\theta)}{\displaystyle\sup_{\theta \in \Theta} L(\theta)}.

Since Θ0⊆Θ\Theta_0 \subseteq \Theta, we have 0≤Λ≤10 \le \Lambda \le 1. Small values of Λ\Lambda indicate that the data are much better explained by some θ\theta outside Θ0\Theta_0 than by any θ\theta inside it, i.e. evidence against H0H_0. The likelihood ratio test at nominal size α\alpha rejects H0H_0 when

Λ≤cα,\Lambda \le c_\alpha,

where cαc_\alpha is chosen (exactly or asymptotically) so that Pr⁡(Λ≤cα∣H0)=α\Pr(\Lambda \le c_\alpha \mid H_0) = \alpha. In the important nested case where Θ\Theta has dimension kk, Θ0\Theta_0 is defined by fixing rr smooth functionally independent constraints (so Θ0\Theta_0 has dimension k−rk-r), the test statistic is usually reported as

G=−2log⁡Λ=2[ℓ(θ^)−ℓ(θ^0)],G = -2\log \Lambda = 2\left[\ell(\hat\theta) - \ell(\hat\theta_0)\right],

where θ^=arg⁡max⁡θ∈Θℓ(θ)\hat\theta = \arg\max_{\theta\in\Theta} \ell(\theta) is the unrestricted MLE and θ^0=arg⁡max⁡θ∈Θ0ℓ(θ)\hat\theta_0 = \arg\max_{\theta \in \Theta_0} \ell(\theta) is the restricted MLE.

Intuition

Think of the log-likelihood function ℓ(θ)\ell(\theta) as a hill over the parameter space. The unrestricted MLE θ^\hat\theta sits at the summit. The null hypothesis Θ0\Theta_0 carves out a lower-dimensional slice of the landscape — a subspace, curve, or hyperplane — and θ^0\hat\theta_0 is the highest point on that slice. If Θ0\Theta_0 happens to pass near the true summit, then ℓ(θ^0)\ell(\hat\theta_0) will be almost as large as ℓ(θ^)\ell(\hat\theta), and the drop in log-likelihood, ℓ(θ^)−ℓ(θ^0)\ell(\hat\theta) - \ell(\hat\theta_0), will be small. If Θ0\Theta_0 cuts through a region far from the summit, the drop is large. The remarkable fact — Wilks's insight — is that near the summit, a smooth log-likelihood behaves like a (multivariate) downward-opening parabola, because the first-order term in a Taylor expansion vanishes at the maximum. The "drop in height" when restricting to a (k−r)(k-r)-dimensional slice through (approximately) the same point is then a quadratic form in Gaussian-distributed estimation errors — exactly the structure that produces a chi-squared distribution with rr degrees of freedom, rr being the codimension of the slice.

Figure
The log-likelihood "hill": unrestricted maximum vs. maximum on the constrained line Θ0\Theta_0
θ̂ (summit) Θ₀ (constraint) θ̂₀ (best on Θ₀) ℓ(θ̂) − ℓ(θ̂₀) ≈ small quadratic gap Level curves of ℓ(θ); gap along Θ₀ drives G = 2[ℓ(θ̂) − ℓ(θ̂₀)]

Central theorem

TheoremWilks's theorem (asymptotic distribution of the generalized likelihood ratio statistic)

Let X1,X2,…X_1, X_2, \dots be i.i.d. with density f(x;θ)f(x;\theta), θ∈Θ⊆Rk\theta \in \Theta \subseteq \mathbb{R}^k an open set. Assume: 1. (Identifiability) θ≠θ′\theta \ne \theta' implies f(⋅;θ)≠f(⋅;θ′)f(\cdot;\theta) \ne f(\cdot;\theta') (as densities/mass functions, a.e.). 2. (Smoothness) The support of f(⋅;θ)f(\cdot;\theta) does not depend on θ\theta, and ℓ1(θ)=log⁡f(X1;θ)\ell_1(\theta) = \log f(X_1;\theta) is three times continuously differentiable in θ\theta in a neighborhood of the true value θ0\theta_0, with third derivatives dominated by an integrable function (regularity conditions sufficient for the standard MLE asymptotics). 3. (Nonsingular Fisher information) The Fisher information matrix

I(θ)=−Eθ ⁣[∂2∂θ ∂θ⊤log⁡f(X1;θ)]I(\theta) = -E_\theta\!\left[\frac{\partial^2}{\partial\theta\,\partial\theta^\top}\log f(X_1;\theta)\right]

exists and is finite and positive definite in a neighborhood of θ0\theta_0. 4. (Interior null) Θ0={θ∈Θ:g(θ)=0}\Theta_0 = \{\theta \in \Theta : g(\theta) = 0\} for a function g:Θ→Rrg:\Theta \to \mathbb{R}^r that is continuously differentiable with Jacobian of full rank rr at θ0\theta_0, and the true parameter θ0\theta_0 lies in Θ0\Theta_0 and is an interior point of Θ\Theta (not on the boundary of Θ\Theta or of Θ0\Theta_0 relative to Θ\Theta). 5. (Consistency) The unrestricted MLE θ^n\hat\theta_n and restricted MLE θ^0,n\hat\theta_{0,n} exist, are consistent for θ0\theta_0, and are interior local maximizers with probability tending to 1. Then, under H0:θ0∈Θ0H_0: \theta_0 \in \Theta_0, as n→∞n \to \infty,

Gn=−2log⁡Λn=2[ℓ(θ^n)−ℓ(θ^0,n)] →d χr2,G_n = -2\log \Lambda_n = 2\left[\ell(\hat\theta_n) - \ell(\hat\theta_{0,n})\right] \ \xrightarrow{d}\ \chi^2_r,

where rr is the number of independent constraints imposed by H0H_0 (the codimension of Θ0\Theta_0 in Θ\Theta).

Proof

We use the standard Taylor-expansion technique for MLE asymptotics, applied twice: once at the unrestricted MLE, once at the restricted MLE, then combined via a change of coordinates that separates the "free" and "constrained" directions.

Step 1: Reduction to a linear constraint via reparametrization. Since gg has Jacobian of rank rr at θ0\theta_0, the implicit function theorem lets us choose (locally, in a neighborhood of θ0\theta_0) a smooth reparametrization ϕ=h(θ)∈Rk\phi = h(\theta) \in \mathbb{R}^k such that hh is a diffeomorphism near θ0\theta_0 and

Θ0 corresponds to {ϕ:ϕk−r+1=⋯=ϕk=0}.\Theta_0 \text{ corresponds to } \{\phi : \phi_{k-r+1} = \cdots = \phi_k = 0\}.

Write ϕ=(ψ,η)\phi = (\psi,\eta) with ψ∈Rk−r\psi \in \mathbb{R}^{k-r} free and η∈Rr\eta \in \mathbb{R}^r the constrained block, so Θ0\Theta_0 becomes {η=0}\{\eta = 0\}. Because likelihood ratios and MLEs are invariant under smooth reparametrization (ℓ\ell evaluated at the corresponding point is unchanged), it suffices to prove the theorem in this new coordinate system with a genuinely linear constraint η=0\eta = 0. Relabel ϕ0=(ψ0,0)\phi_0 = (\psi_0, 0) for the true parameter value in the new coordinates.

Step 2: Asymptotic normality of the unrestricted MLE. Under conditions 2–3, the standard MLE consistency and asymptotic normality theorem (proved via a Taylor expansion of the score equation ∇ℓ(ϕ^n)=0\nabla \ell(\hat\phi_n) = 0 around ϕ0\phi_0) gives

n(ϕ^n−ϕ0)→dN(0, I(ϕ0)−1),\sqrt{n}(\hat\phi_n - \phi_0) \xrightarrow{d} N\big(0,\ I(\phi_0)^{-1}\big),

where I(ϕ0)I(\phi_0) is the Fisher information matrix in the ϕ\phi-coordinates (which exists and is positive definite by condition 3 and the chain rule for Fisher information under diffeomorphisms). This is a prerequisite result from Maximum Likelihood Estimation, not re-derived here, but we use its proof technique — a second-order Taylor expansion of ℓ\ell around the true parameter — again below.

Step 3: Second-order expansion of ℓ\ell around ϕ^n\hat\phi_n. Let ϕ^n=(ψ^n,η^n)\hat\phi_n = (\hat\psi_n,\hat\eta_n) be the unrestricted MLE in the new coordinates. Because ℓ\ell is three times continuously differentiable and ∇ℓ(ϕ^n)=0\nabla \ell(\hat\phi_n) = 0 (first-order condition at an interior maximizer, which holds with probability →1\to 1 by condition 5), Taylor's theorem gives, for any ϕ\phi near ϕ^n\hat\phi_n,

ℓ(ϕ)=ℓ(ϕ^n)−12n(ϕ−ϕ^n)⊤Jn(ϕ−ϕ^n)+Rn(ϕ),\ell(\phi) = \ell(\hat\phi_n) - \tfrac{1}{2} n(\phi - \hat\phi_n)^\top J_n (\phi - \hat\phi_n) + R_n(\phi),

where Jn=−1n∇2ℓ(ϕ^n)J_n = -\tfrac{1}{n}\nabla^2\ell(\hat\phi_n) and Rn(ϕ)=oP(n∥ϕ−ϕ^n∥2)R_n(\phi) = o_P(n\|\phi - \hat\phi_n\|^2) uniformly over ϕ\phi in a shrinking neighborhood of ϕ0\phi_0, by the dominated third-derivative condition (condition 2) together with ϕ^n→ϕ0\hat\phi_n \to \phi_0 in probability (condition 5) and ϕ→ϕ0\phi \to \phi_0. By the law of large numbers and continuity, Jn→pI(ϕ0)J_n \xrightarrow{p} I(\phi_0).

Step 4: Evaluate at the restricted MLE. The restricted MLE ϕ^0,n=(ψ^0,n,0)\hat\phi_{0,n} = (\hat\psi_{0,n}, 0) maximizes ℓ\ell over {η=0}\{\eta = 0\}. Substituting ϕ=(ψ,0)\phi = (\psi,0) into Step 3's expansion and re-maximizing over ψ\psi (equivalently, directly re-expanding around ϕ^n\hat\phi_n restricted to the subspace) yields, after minimizing the quadratic form over the free block ψ\psi,

2[ℓ(ϕ^n)−ℓ(ϕ^0,n)]=n(η^n−0)⊤[(Jn−1)ηη]−1(η^n−0)+oP(1),2\left[\ell(\hat\phi_n) - \ell(\hat\phi_{0,n})\right] = n(\hat\eta_n - 0)^\top \big[ (J_n^{-1})_{\eta\eta} \big]^{-1} (\hat\eta_n - 0) + o_P(1),

where (Jn−1)ηη(J_n^{-1})_{\eta\eta} denotes the r×rr\times r block of Jn−1J_n^{-1} corresponding to the η\eta-coordinates. This is the standard fact that minimizing a positive-definite quadratic form 12v⊤Jv\tfrac12 v^\top J v over the ψ\psi-block, for fixed η\eta-block, produces the Schur-complement quadratic form 12η⊤[(J−1)ηη]−1η\tfrac12 \eta^\top [(J^{-1})_{\eta\eta}]^{-1}\eta in η\eta alone; this is a direct algebraic computation from block matrix inversion, not an appeal to any unproved result.

Step 5: Distribution of the quadratic form. By Step 2 (applied to the η\eta-block, using that marginals of jointly asymptotically normal vectors are asymptotically normal),

n η^n→dN(0, [I(ϕ0)−1]ηη),\sqrt{n}\,\hat\eta_n \xrightarrow{d} N\big(0,\ [I(\phi_0)^{-1}]_{\eta\eta}\big),

since η0=0\eta_0 = 0 under H0H_0. Also Jn→pI(ϕ0)J_n \xrightarrow{p} I(\phi_0) implies (by continuity of matrix inversion and the Schur complement, both continuous maps on the relevant open set of nonsingular matrices)

[(Jn−1)ηη]−1→p[(I(ϕ0)−1)ηη]−1=:Σ−1.\big[(J_n^{-1})_{\eta\eta}\big]^{-1} \xrightarrow{p} \big[(I(\phi_0)^{-1})_{\eta\eta}\big]^{-1} =: \Sigma^{-1}.

Let Zn=n η^nZ_n = \sqrt{n}\,\hat\eta_n. Then

Gn=Zn⊤[(Jn−1)ηη]−1Zn+oP(1)→dZ⊤Σ−1Z,G_n = Z_n^\top \big[(J_n^{-1})_{\eta\eta}\big]^{-1} Z_n + o_P(1) \xrightarrow{d} Z^\top \Sigma^{-1} Z,

by Slutsky's theorem, where Z∼N(0,Σ)Z \sim N(0,\Sigma), applying the continuous mapping theorem to the pair (Zn,[(Jn−1)ηη]−1)(Z_n, [(J_n^{-1})_{\eta\eta}]^{-1}) converging jointly (in distribution for ZnZ_n, in probability for the matrix) to (Z,Σ−1)(Z, \Sigma^{-1}).

Step 6: Identify the quadratic form as χr2\chi^2_r. Since Σ\Sigma is a positive definite r×rr\times r matrix (it is a principal-block inverse of the positive definite matrix I(ϕ0)−1I(\phi_0)^{-1}, hence positive definite), write Σ=Σ1/2Σ1/2\Sigma = \Sigma^{1/2}\Sigma^{1/2} with Σ1/2\Sigma^{1/2} symmetric positive definite. Set W=Σ−1/2ZW = \Sigma^{-1/2}Z. Then W∼N(0,Σ−1/2ΣΣ−1/2)=N(0,Ir)W \sim N(0, \Sigma^{-1/2}\Sigma\Sigma^{-1/2}) = N(0, I_r), and

Z⊤Σ−1Z=W⊤W=∑j=1rWj2,Z^\top \Sigma^{-1} Z = W^\top W = \sum_{j=1}^r W_j^2,

a sum of rr independent squared standard normals — precisely the definition of a χr2\chi^2_r random variable (this is the standard construction of the chi-squared law from Standard Tests: t, Chi-Squared, and F). Combining Steps 5 and 6,

Gn→dχr2,G_n \xrightarrow{d} \chi^2_r,

which completes the proof. ■\blacksquare

∎
CorollaryReduction to classical tests

When f(x;θ)f(x;\theta) is a one-parameter exponential family and H0:θ=θ0H_0: \theta = \theta_0 is a simple hypothesis (r=k=1r = k = 1), the LRT statistic Gn=−2log⁡ΛnG_n = -2\log\Lambda_n is asymptotically equivalent (to first order) to the squared Wald statistic

(θ^n−θ0se^(θ^n))2,\left(\frac{\hat\theta_n - \theta_0}{\widehat{\mathrm{se}}(\hat\theta_n)}\right)^2,

which itself converges to χ12\chi^2_1, matching the square of the classical zz-statistic. This recovers, asymptotically, the familiar one-parameter tests as a special case of Wilks's theorem.

Worked examples

Example 1 — Exponential rate parameter, asymmetric numbers. Suppose lifetimes X1,…,XnX_1,\dots,X_n are i.i.d. Exponential(λ)\mathrm{Exponential}(\lambda) with density f(x;λ)=λe−λxf(x;\lambda) = \lambda e^{-\lambda x}, x>0x>0. We test

H0:λ=λ0=2vs.H1:λ≠2,H_0: \lambda = \lambda_0 = 2 \qquad \text{vs.} \qquad H_1: \lambda \ne 2,

with data n=25n = 25 observations having sample mean Xˉ=0.61\bar X = 0.61 (so ∑Xi=15.25\sum X_i = 15.25). The log-likelihood is ℓ(λ)=nlog⁡λ−λ∑Xi\ell(\lambda) = n\log\lambda - \lambda \sum X_i. The unrestricted MLE solves ℓ′(λ)=n/λ−∑Xi=0\ell'(\lambda) = n/\lambda - \sum X_i = 0, giving

λ^=n∑Xi=2515.25=1.63934…\hat\lambda = \frac{n}{\sum X_i} = \frac{25}{15.25} = 1.63934\ldots

Under H0H_0, Θ0={2}\Theta_0 = \{2\} is a single point, so λ^0=2\hat\lambda_0 = 2 trivially. Here k=1k=1, and Θ0\Theta_0 has dimension 00, so r=1r = 1. Compute the two log-likelihoods (dropping no constants, since they cancel in the difference but we track them for clarity):

ℓ(λ^)=25log⁡(1.63934)−1.63934×15.25=25(0.494697)−25.0000=12.36743−25.00000=−12.63257,\ell(\hat\lambda) = 25\log(1.63934) - 1.63934 \times 15.25 = 25(0.494697) - 25.0000 = 12.36743 - 25.00000 = -12.63257,
ℓ(λ0)=25log⁡(2)−2×15.25=25(0.693147)−30.5=17.32868−30.5=−13.17132.\ell(\lambda_0) = 25\log(2) - 2 \times 15.25 = 25(0.693147) - 30.5 = 17.32868 - 30.5 = -13.17132.

(We used ∑Xiλ^=n\sum X_i \hat\lambda = n exactly, since λ^=n/∑Xi\hat\lambda = n/\sum X_i.) The test statistic is

G=2[ℓ(λ^)−ℓ(λ0)]=2[−12.63257−(−13.17132)]=2(0.53875)=1.0775.G = 2[\ell(\hat\lambda) - \ell(\lambda_0)] = 2[-12.63257 - (-13.17132)] = 2(0.53875) = 1.0775.

Compare to χ12\chi^2_1: the 5%5\% critical value is 3.8413.841. Since 1.0775<3.8411.0775 < 3.841, we fail to reject H0:λ=2H_0:\lambda = 2 at the 5%5\% level; the pp-value is Pr⁡(χ12>1.0775)≈0.2993\Pr(\chi^2_1 > 1.0775) \approx 0.2993, so the observed rate is quite compatible with λ0=2\lambda_0=2. As a sanity check, note the Wald approximation: se^(λ^)=λ^/n=1.63934/5=0.327868\widehat{\mathrm{se}}(\hat\lambda) = \hat\lambda/\sqrt n = 1.63934/5 = 0.327868, and the Wald statistic squared is

(1.63934−20.327868)2=(−0.360660.327868)2=(−1.1002)2=1.2104,\left(\frac{1.63934 - 2}{0.327868}\right)^2 = \left(\frac{-0.36066}{0.327868}\right)^2 = (-1.1002)^2 = 1.2104,

which is close to but not identical to G=1.0775G = 1.0775 — exactly as expected, since the LRT and Wald statistics agree only asymptotically, and n=25n=25 is a moderate sample.

Example 2 — Multinomial goodness-of-fit with unequal cell counts. A genetics-style cross is predicted (under a proposed model H0H_0) to produce four phenotype categories in ratio 9:3:3:19:3:3:1. From n=160n = 160 offspring, the observed counts are

(O1,O2,O3,O4)=(86,35,30,9).(O_1,O_2,O_3,O_4) = (86, 35, 30, 9).

Under the general multinomial model, θ=(p1,p2,p3,p4)\theta = (p_1,p_2,p_3,p_4) with ∑pj=1\sum p_j = 1 is unrestricted (except for the simplex constraint), so k=3k = 3 free parameters. Under H0H_0, pjp_j is fixed at (9/16,3/16,3/16,1/16)(9/16, 3/16, 3/16, 1/16), so Θ0\Theta_0 is a single point: r=3r = 3 constraints. The unrestricted MLE of a multinomial probability is simply the sample proportion, p^j=Oj/n\hat p_j = O_j/n:

p^1=86/160=0.5375,p^2=35/160=0.21875,p^3=30/160=0.1875,p^4=9/160=0.05625.\hat p_1 = 86/160 = 0.5375,\quad \hat p_2 = 35/160 = 0.21875,\quad \hat p_3 = 30/160 = 0.1875,\quad \hat p_4 = 9/160 = 0.05625.

The multinomial log-likelihood (up to the multinomial coefficient, which is the same under both models and cancels) is ℓ(p)=∑jOjlog⁡pj\ell(\mathbf p) = \sum_j O_j \log p_j. So

ℓ(p^)=∑jOjlog⁡(Oj/n).\ell(\hat{\mathbf p}) = \sum_j O_j \log(O_j/n).

Compute term by term (natural log):

86log⁡(0.5375)=86(−0.62098)=−53.4043,86\log(0.5375) = 86(-0.62098) = -53.4043,
35log⁡(0.21875)=35(−1.51983)=−53.1939,35\log(0.21875) = 35(-1.51983) = -53.1939,
30log⁡(0.1875)=30(−1.67398)=−50.2194,30\log(0.1875) = 30(-1.67398) = -50.2194,
9log⁡(0.05625)=9(−2.87795)=−25.9016.9\log(0.05625) = 9(-2.87795) = -25.9016.

Sum: ℓ(p^)=−53.4043−53.1939−50.2194−25.9016=−182.7192\ell(\hat{\mathbf p}) = -53.4043 -53.1939 -50.2194 -25.9016 = -182.7192. Under H0H_0, pj0=(9/16,3/16,3/16,1/16)=(0.5625,0.1875,0.1875,0.0625)p_j^0 = (9/16, 3/16, 3/16, 1/16) = (0.5625, 0.1875, 0.1875, 0.0625):

86log⁡(0.5625)=86(−0.57536)=−49.4810,86\log(0.5625) = 86(-0.57536) = -49.4810,
35log⁡(0.1875)=35(−1.67398)=−58.5893,35\log(0.1875) = 35(-1.67398) = -58.5893,
30log⁡(0.1875)=30(−1.67398)=−50.2194,30\log(0.1875) = 30(-1.67398) = -50.2194,
9log⁡(0.0625)=9(−2.77259)=−24.9533.9\log(0.0625) = 9(-2.77259) = -24.9533.

Sum: ℓ(p0)=−49.4810−58.5893−50.2194−24.9533=−183.2430\ell(\mathbf p^0) = -49.4810 - 58.5893 - 50.2194 - 24.9533 = -183.2430. The LRT statistic is

G=2[ℓ(p^)−ℓ(p0)]=2[−182.7192−(−183.2430)]=2(0.5238)=1.0476.G = 2[\ell(\hat{\mathbf p}) - \ell(\mathbf p^0)] = 2[-182.7192 - (-183.2430)] = 2(0.5238) = 1.0476.

Comparing to χ32\chi^2_3 (critical value at 5%5\% is 7.8157.815), G=1.0476G = 1.0476 gives p-value=Pr⁡(χ32>1.0476)≈0.7896p\text{-value} = \Pr(\chi^2_3 > 1.0476) \approx 0.7896; the 9:3:3:19:3:3:1 model fits comfortably. This matches the classical Pearson chi-squared statistic

X2=∑j(Oj−Ej)2Ej,Ej=npj0=(90,30,30,10),X^2 = \sum_j \frac{(O_j - E_j)^2}{E_j}, \qquad E_j = n p_j^0 = (90, 30, 30, 10),

which computes to

X2=(86−90)290+(35−30)230+(30−30)230+(9−10)210=1690+2530+0+110=0.1778+0.8333+0+0.1=1.1111,X^2 = \frac{(86-90)^2}{90} + \frac{(35-30)^2}{30} + \frac{(30-30)^2}{30} + \frac{(9-10)^2}{10} = \frac{16}{90} + \frac{25}{30} + 0 + \frac{1}{10} = 0.1778 + 0.8333 + 0 + 0.1 = 1.1111,

close to G=1.0476G = 1.0476 — the two statistics are asymptotically equivalent (both converge to χ32\chi^2_3 under H0H_0), a second-order Taylor-expansion fact connecting the LRT statistic to Pearson's statistic in multinomial models.

CounterexampleWhy the "interior null" hypothesis is essential

Wilks's theorem requires that the true parameter θ0\theta_0 be an interior point of Θ0\Theta_0's ambient space relative to Θ\Theta — i.e., the null hypothesis must not force the parameter to the boundary of the parameter space. Drop this hypothesis and the chi-squared limit fails.

Concrete failure. Let X1,…,XnX_1,\dots,X_n be i.i.d. N(μ,1)N(\mu,1) with known variance 11, and suppose the parameter space is restricted on physical grounds to Θ=[0,∞)\Theta = [0,\infty) (e.g. μ\mu represents a concentration or an intensity that cannot be negative). Test

H0:μ=0vs.H1:μ>0,H_0: \mu = 0 \qquad \text{vs.}\qquad H_1: \mu > 0,

with μ=0\mu = 0 sitting on the boundary of Θ\Theta. The unrestricted MLE over Θ=[0,∞)\Theta = [0,\infty) is μ^n=max⁡(Xˉn,0)\hat\mu_n = \max(\bar X_n, 0), since the unconstrained maximizer of the Gaussian log-likelihood is Xˉn\bar X_n, projected onto [0,∞)[0,\infty) when Xˉn<0\bar X_n < 0. The restricted MLE under H0H_0 is simply μ^0,n=0\hat\mu_{0,n} = 0. The likelihood ratio statistic is

Gn=2[ℓ(μ^n)−ℓ(0)]={nXˉn2,Xˉn≥0,0,Xˉn<0.G_n = 2\left[\ell(\hat\mu_n) - \ell(0)\right] = \begin{cases} n\bar X_n^2, & \bar X_n \ge 0,\\[4pt] 0, & \bar X_n < 0. \end{cases}

Under H0H_0 (μ0=0\mu_0 = 0), n Xˉn→dN(0,1)\sqrt n\,\bar X_n \xrightarrow{d} N(0,1). So with probability →1/2\to 1/2 (whenever Xˉn<0\bar X_n < 0), GnG_n is exactly zero, and with probability →1/2\to 1/2, Gn=(nXˉn)2G_n = (\sqrt n \bar X_n)^2 converges to the square of a standard normal, i.e. χ12\chi^2_1. The limiting distribution of GnG_n is therefore

Gn→d12δ0+12χ12,G_n \xrightarrow{d} \tfrac12 \delta_0 + \tfrac12 \chi^2_1,

a 50–50 mixture of a point mass at zero and a χ12\chi^2_1 — not χ12\chi^2_1 itself. Using the naive χ12\chi^2_1 critical value 3.8413.841 for a nominal 5%5\% test actually yields an asymptotic size of only 12Pr⁡(χ12>3.841)=0.025\tfrac12 \Pr(\chi^2_1 > 3.841) = 0.025, half the intended level. This single example (due essentially to Chernoff, 1954) shows precisely what the interior-null hypothesis buys you: without it, the quadratic Taylor expansion around θ^0\hat\theta_0 is no longer valid on both sides, because the restricted maximizer is forced to the edge and cannot move freely in all directions — the "free directions" argument in Step 4 of the proof breaks down, since ψ\psi (here, the sign-constrained μ\mu) cannot range over an open neighborhood of ϕ0\phi_0.

Scope, edge cases & common traps

SituationWhat happensFix / correct approach
Null parameter on boundary of Θ\ThetaLimit is a mixture of chi-squareds, not pure χr2\chi^2_rUse Chernoff/Self-Korn theory for boundary problems, or reparametrize
Non-identifiable model under H0H_0 (e.g. testing number of mixture components)Fisher information degenerates; limit is a supremum of a Gaussian process, not χ2\chi^2Use specialized theory (e.g. Davies' test) for mixture order selection
Small nnχr2\chi^2_r is only an asymptotic approximation; exact small-sample distribution may differ substantiallyUse exact tests (e.g. exact Poisson/binomial LRTs) or bootstrap the null distribution
Nuisance parameters present in both modelsStill fine — rr counts only the constrained dimensions, nuisance parameters are "profiled out" automaticallyCompute θ^\hat\theta, θ^0\hat\theta_0 by maximizing over all free coordinates including nuisance ones
Non-nested modelsΘ0⊈Θ\Theta_0 \not\subseteq \Theta — the ratio is not even guaranteed to lie in [0,1][0,1]Use Vuong's test or information criteria (AIC/BIC) instead

The mathematician behind it

Origins · 1 minSamuel S. Wilks and how this was discovered
SS
Samuel S. Wilks
1906–1964 · Princeton University, mathematical statistics

Samuel Stanley Wilks was an American mathematical statistician whose 1938 paper "The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses" established the theorem that bears his name. Wilks worked at Princeton for his entire career, where he built one of the first serious American graduate programs in mathematical statistics, trained influential students (including John Tukey), and was a founding force behind the Annals of Mathematical Statistics. His broader body of work — on multivariate analysis, statistical tables, and the mathematical foundations of experimental design — helped establish statistics as a rigorous mathematical discipline in the United States in the mid-20th century, moving it beyond its earlier home in applied biology and agriculture.

History
- 1928 — Jerzy Neyman and Egon Pearson begin developing hypothesis-testing theory that will culminate in the Neyman-Pearson Lemma.
- 1933 — Neyman and Pearson publish the fundamental lemma establishing the likelihood ratio's optimality for simple-vs-simple tests.
- 1938 — Samuel Wilks proves the asymptotic χr2\chi^2_r distribution of the generalized likelihood ratio statistic for composite hypotheses.
- 1954 — Herman Chernoff analyzes the boundary case, showing the mixture-of-chi-squared phenomenon when the null lies on the boundary of the parameter space.
- 1943–1946 — Abraham Wald develops the closely related Wald test, and C.R. Rao later develops the score (Lagrange multiplier) test, completing the "holy trinity" of asymptotically equivalent tests built on the likelihood.
Why is −2log⁡Λ-2\log\Lambda used instead of Λ\Lambda itself?

Taking −2log⁡Λ-2\log\Lambda transforms a ratio bounded in [0,1][0,1] (awkward to calibrate directly) into a nonnegative quantity that, by Wilks's theorem, has a clean, sample-size-independent asymptotic distribution: χr2\chi^2_r. The factor of 22 is not arbitrary — it exactly cancels the 12\tfrac12 that appears from the second-order Taylor expansion of the log-likelihood around its maximum (where the Hessian contributes a factor of 12\tfrac12 to the quadratic approximation). Without the 22, the limiting distribution would be a chi-squared variable scaled by 12\tfrac12, which is a needless complication; with it, tables and software can report the statistic directly against standard χ2\chi^2 quantiles.

How does the likelihood ratio test relate to the Neyman-Pearson lemma?

The Neyman-Pearson Lemma proves that for testing one simple hypothesis θ=θ0\theta=\theta_0 against another simple alternative θ=θ1\theta=\theta_1, the likelihood ratio f(x;θ0)/f(x;θ1)f(x;\theta_0)/f(x;\theta_1) gives the uniformly most powerful test at any level α\alpha. The generalized likelihood ratio test extends this idea to composite hypotheses, where Θ0\Theta_0 and Θ∖Θ0\Theta\setminus\Theta_0 are not single points but sets, by replacing each hypothesis's likelihood with its supremum over the relevant parameter set. The GLRT is generally not exactly optimal (no uniformly most powerful test typically exists for composite alternatives), but it retains excellent large-sample power properties and reduces to the exact Neyman-Pearson test in the simple-vs-simple special case.

Is the likelihood ratio test the same as the Wald test or the score test?

No, but they are asymptotically equivalent under the null hypothesis and under Wilks's regularity conditions. The Wald test uses only the unrestricted MLE and its estimated standard error (a quadratic approximation around θ^\hat\theta); the score (Lagrange multiplier) test uses only the restricted MLE and the score function evaluated there (a quadratic approximation around θ^0\hat\theta_0); the LRT uses both fitted models. All three converge to the same χr2\chi^2_r limit, but they can differ noticeably in finite samples, as seen in Example 1 above where the Wald statistic (1.211.21) and LRT statistic (1.081.08) were close but not identical.

Does the likelihood ratio test require nested models?

Yes, in its classical form: Θ0\Theta_0 must be a subset of Θ\Theta for the ratio Λ\Lambda to be well-defined and bounded by 11, and Wilks's theorem is stated and proved specifically for this nested setting. For comparing non-nested models (e.g., a lognormal model versus a gamma model for the same data), the ratio of maximized likelihoods is still computable, but it is no longer bounded by 11 and Wilks's χ2\chi^2 limit does not apply. Non-nested comparisons instead use tools like Vuong's test, the Akaike Information Criterion, or the Bayesian Information Criterion, which are built on related but distinct theoretical foundations.

What role does the number of degrees of freedom rr play in practice, and how do I count it correctly?

rr is the number of independent scalar equality constraints that H0H_0 imposes relative to the unrestricted model — equivalently, dim⁡(Θ)−dim⁡(Θ0)\dim(\Theta) - \dim(\Theta_0). In a regression setting, testing whether 33 specific coefficients are simultaneously zero gives r=3r=3; testing a single coefficient gives r=1r=1. In the multinomial example, testing that all four cell probabilities equal specified values under a model with 33 free probabilities (since they must sum to 11) gives r=3r=3, not 44. Miscounting rr — a very common error — leads to comparing GG against the wrong chi-squared distribution and hence to an incorrect pp-value.

What happens to the likelihood ratio test in exact, non-asymptotic Gaussian linear models?

In the classical linear model Y=Xβ+εY = X\beta + \varepsilon with ε∼N(0,σ2I)\varepsilon \sim N(0,\sigma^2 I), testing a linear restriction on β\beta (with σ2\sigma^2 unknown and estimated) produces a likelihood ratio statistic that is a strictly increasing function of the classical FF-statistic from Standard Tests: t, Chi-Squared, and F. In this special case the exact finite-sample null distribution of a monotone transform of Λ\Lambda is known exactly (an Fr, n−kF_{r,\,n-k} distribution), rather than only asymptotically χr2\chi^2_r — a rare instance where the LRT machinery yields an exact, not merely asymptotic, classical test.

Can the likelihood ratio test statistic ever be negative, or can GG come out negative in practice?

No: since Θ0⊆Θ\Theta_0 \subseteq \Theta, the unrestricted maximum likelihood can never be smaller than the restricted maximum likelihood, so Λ≤1\Lambda \le 1 and G=−2log⁡Λ≥0G = -2\log\Lambda \ge 0 always, by construction. If a computed GG comes out negative in an implementation, this signals a numerical or optimization error — most commonly, that the "unrestricted" optimizer failed to find the true global maximum (got stuck in a local optimum with lower likelihood than the constrained fit), which is a common practical pitfall with iterative MLE routines in nonlinear or high-dimensional models.

Before this

Where this leads