Goodness-of-Fit & Empirical Distribution Tests: Chi-Square, Kolmogorov-Smirnov & Lilliefors

Let $X_1, X_2, \dots, X_n$ be i.i.d. observations from a distribution $F$, and suppose we partition the sample space $\mathcal{X}$ into $k$ mutually exclusive and exhaustive cells (categories or intervals) $C_1, C_2, \dots, C_k$ such that $\bigcup_{j=1}^k C_j = \mathcal{X}$ and $C_i \cap C_j = \emptyset$ for $i \ne j$. Let $O_j = \sum_{i=1}^n \mathbb{I}(X_i \in C_j)$ denote the observed frequency of sample observations falling into cell $C_j$, so that $\sum_{j=1}^k O_j = n$. Under a simple null hypothesis $H_0 : F = F_0$, each observation falls into cell $C_j$ with known theoretical probability:

Consequently, under $H_0$, the vector of observed cell counts $\mathbf{O} = (O_1, O_2, \dots, O_k)^\top$ follows an exact Multinomial distribution $\text{Multinomial}(n; p_{10}, p_{20}, \dots, p_{k0})$ with expected cell frequencies $E_j = \mathbb{E}_{H_0}[O_j] = n p_{j0}$. In 1900, Karl Pearson proposed measuring the discrepancy between the observed counts $O_j$ and expected counts $E_j$ via the weighted sum of squared residuals:

Why does division by $E_j = n p_{j0}$ convert a singular multinomial covariance matrix into an exact asymptotic $\chi^2_{k-1}$ distribution? Let us prove Pearson's theorem using multivariate normal asymptotics and projection matrix algebra.

Theorem 6.1 (Pearson's Asymptotic Chi-Square Goodness-of-Fit Theorem)

Under the simple null hypothesis $H_0 : p_j = p_{j0} > 0$ ($j = 1, \dots, k$) with $\sum_{j=1}^k p_{j0} = 1$, as the sample size $n \to \infty$, Pearson's statistic $X^2$ converges in distribution to a central Chi-Square distribution with $k - 1$ degrees of freedom:

Furthermore, the Generalized Likelihood Ratio Test (GLRT) statistic $G^2 = -2 \ln \Lambda = 2 \sum_{j=1}^k O_j \ln\left(\frac{O_j}{E_j}\right)$ is asymptotically equivalent to $X^2$ under $H_0$: $G^2 - X^2 \xrightarrow{P} 0$, so $G^2 \xrightarrow{d} \chi^2_{k-1}$.

Proof:

For each observation $i \in \{1, \dots, n\}$, define the $k$-dimensional standardized indicator vector $\mathbf{Y}_i = (Y_{i1}, \dots, Y_{ik})^\top$, where $Y_{ij} = \frac{\mathbb{I}(X_i \in C_j) - p_{j0}}{\sqrt{p_{j0}}}$. Under $H_0$, $\mathbb{E}[\mathbf{Y}_i] = \mathbf{0}$, and the covariance matrix $\boldsymbol{\Sigma} = \text{Cov}(\mathbf{Y}_i)$ has diagonal entries $\Sigma_{jj} = \frac{p_{j0}(1 - p_{j0})}{p_{j0}} = 1 - p_{j0}$ and off-diagonal entries $\Sigma_{jl} = \frac{-p_{j0} p_{l0}}{\sqrt{p_{j0} p_{l0}}} = -\sqrt{p_{j0} p_{l0}}$ for $j \ne l$. In matrix notation, define the unit vector $\mathbf{u} = (\sqrt{p_{10}}, \sqrt{p_{20}}, \dots, \sqrt{p_{k0}})^\top \in \mathbb{R}^k$, which has Euclidean norm $\|\mathbf{u}\|^2 = \mathbf{u}^\top \mathbf{u} = \sum_{j=1}^k p_{j0} = 1$. Then the covariance matrix $\boldsymbol{\Sigma}$ takes the clean rank-$(k-1)$ projection form:

Notice that $\boldsymbol{\Sigma}$ is symmetric and idempotent:

Now consider the sample sum vector $\mathbf{Z}_n = \frac{1}{\sqrt{n}} \sum_{i=1}^n \mathbf{Y}_i$, whose $j$-th coordinate is precisely $Z_{nj} = \frac{O_j - n p_{j0}}{\sqrt{n p_{j0}}}$. Notice that Pearson's statistic is simply the squared Euclidean norm of $\mathbf{Z}_n$:

By the Multivariate Central Limit Theorem, as $n \to \infty$, $\mathbf{Z}_n \xrightarrow{d} \mathbf{Z} \sim \mathcal{N}_k(\mathbf{0}, \boldsymbol{\Sigma})$, where $\boldsymbol{\Sigma} = \mathbf{I}_k - \mathbf{u}\mathbf{u}^\top$ is an idempotent projection matrix of rank equal to its trace:

By the spectral theorem for quadratic forms in normal variables (since $\mathbf{Z}^\top \mathbf{Z} = \mathbf{Z}^\top \boldsymbol{\Sigma} \mathbf{Z}$ with probability $1$ because $\mathbf{u}^\top \mathbf{Z}_n = \frac{1}{\sqrt{n}}\sum_{j=1}^k (O_j - n p_{j0}) = 0$), $\mathbf{Z}_n^\top \mathbf{Z}_n \xrightarrow{d} \mathbf{Z}^\top \mathbf{Z} \sim \chi^2_{\text{rank}(\boldsymbol{\Sigma})} = \chi^2_{k-1}$. Finally, expanding $G^2 = 2 \sum_{j=1}^k O_j \ln\left(1 + \frac{O_j - E_j}{E_j}\right)$ via a second-order Taylor series $\ln(1 + \delta) = \delta - \frac{1}{2}\delta^2 + O(\delta^3)$ gives $G^2 = 2 \sum_{j=1}^k (E_j + (O_j - E_j))\left[\frac{O_j - E_j}{E_j} - \frac{(O_j - E_j)^2}{2 E_j^2} + O_p(n^{-3/2})\right] = \sum_{j=1}^k \frac{(O_j - E_j)^2}{E_j} + o_p(1) = X^2 + o_p(1)$, completing the proof.

What happens when the null hypothesis $H_0 : p_j = p_j(\boldsymbol{\theta})$ is composite, depending on an unknown $s$-dimensional parameter vector $\boldsymbol{\theta} = (\theta_1, \dots, \theta_s)^\top \in \Theta \subseteq \mathbb{R}^s$ ($s < k - 1$)? Karl Pearson originally believed in 1900 that replacing $p_{j0}$ by $p_j(\hat{\boldsymbol{\theta}})$ did not alter the asymptotic $\chi^2_{k-1}$ distribution! R. A. Fisher (1922, 1924) proved Pearson wrong in a monumental clash of statistical titans, showing that estimating $s$ independent parameters by multinomial Maximum Likelihood (or Minimum Chi-Square) projects the standardized residual vector $\hat{\mathbf{Z}}_n$ onto a lower-dimensional subspace of dimension $k - 1 - s$.

Theorem 6.2 (Fisher's Degrees-of-Freedom Reduction Theorem and Chernoff-Lehmann Phenomenon)

Suppose the cell probabilities $p_j(\boldsymbol{\theta}) > 0$ ($j = 1, \dots, k$) are continuously differentiable functions of $\boldsymbol{\theta} \in \Theta \subseteq \mathbb{R}^s$ ($s < k - 1$) such that the $k \times s$ Jacobian matrix $\mathbf{B}(\boldsymbol{\theta})$ with entries $B_{jr}(\boldsymbol{\theta}) = \frac{1}{\sqrt{p_j(\boldsymbol{\theta})}} \frac{\partial p_j(\boldsymbol{\theta})}{\partial \theta_r}$ has full column rank $s$.

1. If $\hat{\boldsymbol{\theta}}_{\text{MLE}}$ is the multinomial (grouped-data) Maximum Likelihood Estimator maximizing $\sum_{j=1}^k O_j \ln p_j(\boldsymbol{\theta})$ (or the Minimum Chi-Square estimator minimizing $\sum_{j=1}^k \frac{(O_j - n p_j(\boldsymbol{\theta}))^2}{n p_j(\boldsymbol{\theta})}$), then as $n \to \infty$:

1. (Chernoff-Lehmann Theorem, 1954): If instead $\tilde{\boldsymbol{\theta}}_n$ is estimated from the ungrouped raw data $X_1, \dots, X_n$ before binning, then $X^2(\tilde{\boldsymbol{\theta}}_n)$ converges in distribution to a weighted sum of independent $\chi^2_1$ variables:

Hence the limiting distribution of $X^2(\tilde{\boldsymbol{\theta}}_n)$ is stochastically bounded between $\chi^2_{k - 1 - s}$ and $\chi^2_{k - 1}$.

In practice, to ensure the validity of the multinomial normal approximation, Cochran's Rule (1952, 1954) recommends choosing bins so that every expected cell frequency is at least $1$ ($E_j \ge 1$) and at least $80\%$ of the cells satisfy $E_j \ge 5$; if tail cells violate this criterion, adjacent tail cells should be pooled (merging cells reduces $k$ accordingly).

Pearson's $\chi^2$ framework extends naturally to bivariate and multivariate categorical data arranged in contingency tables. Suppose a random sample of size $n$ is cross-classified according to two categorical attributes: Factor $A$ with $R$ levels (rows $i = 1, \dots, R$) and Factor $B$ with $C$ levels (columns $j = 1, \dots, C$), forming an $R \times C$ contingency table with observed cell counts $O_{ij}$, row totals $R_i = \sum_{j=1}^C O_{ij}$, column totals $C_j = \sum_{i=1}^R O_{ij}$, and grand total $n = \sum_{i=1}^R \sum_{j=1}^C O_{ij}$.

Two distinct experimental sampling designs lead to the exact same algebraic test statistic:

1. Test of Independence (Single Multinomial Sample of size $n$): A single random sample of size $n$ is drawn from a bivariate population where $\pi_{ij} = \mathbb{P}(\text{Row } i, \text{Col } j)$, $\pi_{i\cdot} = \sum_{j=1}^C \pi_{ij}$, and $\pi_{\cdot j} = \sum_{i=1}^R \pi_{ij}$. We test the null hypothesis of statistical independence:

Here there are $k = RC$ cells, so the unconstrained multinomial has $RC - 1$ degrees of freedom. Under $H_0$, we estimate $R - 1$ independent row marginals ($\hat{\pi}_{i\cdot} = R_i / n$) and $C - 1$ independent column marginals ($\hat{\pi}_{\cdot j} = C_j / n$), so $s = (R - 1) + (C - 1)$ parameters are estimated. By Fisher's theorem, the degrees of freedom are:

2. Test of Homogeneity ($R$ Independent Multinomial Samples of Fixed Row Sizes $R_1, \dots, R_R$): Independent samples of fixed sizes $R_i$ are drawn from $R$ distinct subpopulations, and we test whether the conditional column distributions are identical across all subpopulations:

Under both designs, the MLE of the expected cell count under $H_0$ is $\hat{E}_{ij} = \frac{R_i C_j}{n}$, and Pearson's $\chi^2$ statistic is:

In the special case of a $2 \times 2$ contingency table with cell counts $(O_{11}, O_{12}; O_{21}, O_{22})$, Pearson's statistic simplifies algebraically to:

To improve the continuous $\chi^2_1$ approximation to the discrete hypergeometric conditional distribution (Fisher's exact test) when sample sizes are moderate, Frank Yates (1934) introduced Yates' Continuity Correction, which subtracts $n/2$ from the absolute cross-product difference $|O_{11} O_{22} - O_{12} O_{21}|$ (equivalent to subtracting $0.5$ from each $|O_{ij} - \hat{E}_{ij}|$):

When observational studies are confounded by a third stratifying variable $Z$ with $K$ strata ($k = 1, \dots, K$, such as hospital centers or age brackets), pooling the $K$ tables crudely into a single marginal $2 \times 2$ table can induce Simpson's Paradox. To test conditional independence $H_0 : \psi_1 = \psi_2 = \dots = \psi_K = 1$ across $K$ stratified $2 \times 2$ tables (where $\psi_k = \frac{\pi_{11k}\pi_{22k}}{\pi_{12k}\pi_{21k}}$ is the odds ratio in stratum $k$), William Cochran (1954) and Nathan Mantel and William Haenszel (1959) developed the Cochran-Mantel-Haenszel (CMH) Test. Conditional on the margins $(R_{1k}, R_{2k}, C_{1k}, C_{2k}, n_k)$ in stratum $k$, the upper-left cell count $O_{11k}$ follows a central Hypergeometric distribution under $H_0$ with exact mean and variance:

Summing the deviations $O_{11k} - \mu_{11k}$ across all $K$ independent strata yields the Cochran-Mantel-Haenszel statistic:

where $c_{\text{cc}} = 0.5$ (with continuity correction) or $0$ (uncorrected), alongside the Mantel-Haenszel Common Odds Ratio Estimator $\hat{\psi}_{\text{MH}} = \frac{\sum_{k=1}^K O_{11k} O_{22k} / n_k}{\sum_{k=1}^K O_{12k} O_{21k} / n_k}$.

While Pearson's $\chi^2$ test is versatile for categorical and discrete data, applying it to continuous random variables $X_1, \dots, X_n$ suffers from two serious theoretical drawbacks: (1) arbitrary choice of bin boundaries $C_1, \dots, C_k$ can alter the test conclusion, and (2) grouping continuous observations into coarse intervals discards exact intra-bin ordering information and destroys statistical power against local alternatives. To preserve full information in continuous samples, we construct the Empirical Distribution Function (EDF) $F_n(x)$, defined for every real number $x \in \mathbb{R}$ as the proportion of sample observations less than or equal to $x$:

In terms of the order statistics $X_{(1)} \le X_{(2)} \le \dots \le X_{(n)}$, $F_n(x)$ is a right-continuous step function (a càdlàg function in the Skorokhod space $D[-\infty, \infty]$) that jumps upward by $1/n$ at each distinct observation:

For any fixed real number $x \in \mathbb{R}$, the indicator variable $\mathbb{I}(X_i \le x) \stackrel{\text{i.i.d.}}{\sim} \text{Bernoulli}(F(x))$, so the count $n F_n(x) \sim \text{Binomial}(n, F(x))$. Consequently, $F_n(x)$ is an unbiased estimator of $F(x)$ with exact pointwise moments and covariance structure:

• Pointwise Mean and Variance: $\mathbb{E}[F_n(x)] = F(x)$, and $\text{Var}(F_n(x)) = \frac{F(x)[1 - F(x)]}{n}$.
• Exact Covariance Function (for any $x, y \in \mathbb{R}$):

• Pointwise Strong Consistency and Asymptotic Normality: By Kolmogorov's Strong Law of Large Numbers, $F_n(x) \xrightarrow{\text{a.s.}} F(x)$ for each fixed $x$, and by the Lindeberg-Lévy Central Limit Theorem, $\sqrt{n}(F_n(x) - F(x)) \xrightarrow{d} \mathcal{N}\left(0, F(x)[1 - F(x)]\right)$.

Remarkably, Valerio Glivenko and Francesco Paolo Cantelli (1933) proved that the pointwise almost-sure convergence of $F_n(x)$ to $F(x)$ is automatically uniform across the entire real line $\mathbb{R}$, a result celebrationally termed the Fundamental Theorem of Mathematical Statistics.

Theorem 6.3 (The Glivenko-Cantelli Theorem)

Let $X_1, X_2, \dots, X_n$ be i.i.d. random variables with common cumulative distribution function $F(x)$. Define the Kolmogorov supremum distance $D_n = \|F_n - F\|_\infty = \sup_{x \in \mathbb{R}} |F_n(x) - F(x)|$. Then as $n \to \infty$, $F_n$ converges uniformly to $F$ almost surely:

Proof:

Fix any integer $m \ge 1$, and partition the probability interval $[0, 1]$ into $m$ subintervals of length $1/m$ via the generalized quantile points $-\infty = x_0 < x_1 \le x_2 \le \dots \le x_{m-1} < x_m = +\infty$, where $x_j = \inf\{x \in \mathbb{R} : F(x) \ge j/m\}$ for $j = 1, \dots, m-1$. By right-continuity and monotonicity of $F$, on each open interval $(x_{j-1}, x_j)$, the oscillation of $F$ cannot exceed $1/m$: $F(x_j^-) - F(x_{j-1}) \le \frac{1}{m}$. Now let $x \in \mathbb{R}$ be arbitrary; if $x$ lies in the open interval $(x_{j-1}, x_j)$ for some $j \in \{1, \dots, m\}$, monotonicity of both $F_n$ and $F$ implies:

Since $F(x_j^-) \le F(x_{j-1}) + \frac{1}{m}$, both bounds are controlled by the finite grid of endpoints $\{x_1, \dots, x_{m-1}\}$:

Taking the supremum over all $x \in \mathbb{R}$ yields the deterministic sandwich inequality:

Because $\{x_1, \dots, x_{m-1}\}$ is a finite set of $2(m-1)$ evaluation points, applying the classical Strong Law of Large Numbers to each of the $2(m-1)$ points guarantees that with probability one, $\lim_{n \to \infty} \max_{1 \le j \le m-1} \max\{|F_n(x_j) - F(x_j)|, |F_n(x_j^-) - F(x_j^-)|\} = 0$. Therefore, with probability one:

Since $m \ge 1$ was completely arbitrary, taking a countable intersection over $m = 1, 2, \dots$ yields $\mathbb{P}(\lim_{n \to \infty} D_n = 0) = 1$.

Moreover, Monroe Donsker (1952) proved the Functional Central Limit Theorem for the Empirical Process: when $F$ is continuous, the scaled empirical process $\alpha_n(x) = \sqrt{n}\left(F_n(x) - F(x)\right)$, transformed via $t = F(x) \in [0, 1]$ into the uniform empirical process $U_n(t) = \sqrt{n}\left(G_n(t) - t\right)$ on $[0, 1]$, converges weakly in $D[0, 1]$ to a Standard Brownian Bridge (pinned Wiener process) $\{\mathbb{B}(t) : 0 \le t \le 1\}$:

where $\{W(t) : t \ge 0\}$ is a standard Brownian motion, and $\mathbb{B}(t)$ is a continuous zero-mean Gaussian process with $\mathbb{B}(0) = \mathbb{B}(1) = 0$ and covariance kernel $\text{Cov}(\mathbb{B}(s), \mathbb{B}(t)) = \min(s, t) - st$. Even in finite samples, the Massart (1990) sharp Dvoretzky-Kiefer-Wolfowitz (DKW) Inequality provides an exact non-asymptotic concentration bound for all $n \ge 1$ and $\epsilon > 0$:

Armed with the Empirical Distribution Function $F_n(x)$, Andrey Kolmogorov (1933) and Nikolai Smirnov (1939) constructed exact distribution-free goodness-of-fit tests based on the $L^\infty$ supremum norm between $F_n(x)$ and a continuous null CDF $F_0(x)$. For testing the simple hypothesis $H_0 : F(x) = F_0(x)$ for all $x \in \mathbb{R}$ (where $F_0$ is completely specified and continuous), define the one-sided and two-sided Kolmogorov-Smirnov (K-S) statistics:

Notice that $D_n^+$ tests against the one-sided stochastic dominance alternative $H_1^+ : F(x) \ge F_0(x)$ (values of $X$ stochastically smaller than $F_0$), whereas $D_n^-$ tests against $H_1^- : F(x) \le F_0(x)$ (values of $X$ stochastically larger than $F_0$). How do we compute $D_n^+$, $D_n^-$, and $D_n$ from a finite sample $X_1, \dots, X_n$ without taking an uncountable supremum over $x \in \mathbb{R}$? Because $F_n(x)$ is constant equal to $(i-1)/n$ on each open interval $(X_{(i-1)}, X_{(i)})$ while $F_0(x)$ is monotonically non-decreasing, the differences $F_n(x) - F_0(x)$ and $F_0(x) - F_n(x)$ achieve their extrema strictly at the order statistics $X_{(1)} \le X_{(2)} \le \dots \le X_{(n)}$:

Theorem 6.4 (Distribution-Free Property and Kolmogorov-Smirnov Limiting Distributions)

Let $X_1, \dots, X_n$ be i.i.d. from a continuous CDF $F_0(x)$.

1. Exact Finite-Sample Distribution-Freeness: The exact distributions of $D_n^+$, $D_n^-$, and $D_n$ under $H_0$ are completely independent of $F_0$, and equal the distributions of the uniform K-S statistics for $U_1, \dots, U_n \stackrel{\text{i.i.d.}}{\sim} \text{Uniform}(0, 1)$.

2. Smirnov's Asymptotic Distribution for One-Sided $D_n^+$ and $D_n^-$: For every $\lambda > 0$:

3. Kolmogorov's Asymptotic Distribution for Two-Sided $D_n$: For every $\lambda > 0$, $\sqrt{n} D_n \xrightarrow{d} \sup_{0 \le t \le 1} |\mathbb{B}(t)|$, whose exact CDF is the alternating theta-function series:

Proof of Distribution-Freeness and Asymptotic Convergence:

By the Probability Integral Transform, when $X_i \sim F_0$ with $F_0$ continuous, the transformed random variables $U_i = F_0(X_i)$ are i.i.d. $\text{Uniform}(0, 1)$, and their order statistics are $U_{(i)} = F_0(X_{(i)})$. Substituting $U_{(i)} = F_0(X_{(i)})$ into the computational formulas above yields:

Because the joint distribution of $(U_{(1)}, \dots, U_{(n)})$ is the Dirichlet/uniform order-statistic distribution on the simplex $0 < u_1 < \dots < u_n < 1$ with density $n!$ independent of $F_0$, $D_n^+$, $D_n^-$, and $D_n$ are strictly distribution-free! Furthermore, by Donsker's Theorem and the Continuous Mapping Theorem applied to the continuous functionals $g^+(f) = \sup_{t \in [0, 1]} f(t)$ and $g(f) = \sup_{t \in [0, 1]} |f(t)|$ on $D[0, 1]$, we have $\sqrt{n} D_n^+ \xrightarrow{d} \sup_{0 \le t \le 1} \mathbb{B}(t)$ and $\sqrt{n} D_n \xrightarrow{d} \sup_{0 \le t \le 1} |\mathbb{B}(t)|$. Applying the reflection principle of André to the pinned Brownian bridge $\mathbb{B}(t)$ across the single barrier $+\lambda$ yields $\mathbb{P}(\sup_{0 \le t \le 1} \mathbb{B}(t) \ge \lambda) = e^{-2\lambda^2}$, while repeated inclusion-exclusion reflections across the double barriers $\pm \lambda$ yield Kolmogorov's alternating series $1 - K(\lambda) = 2 \sum_{j=1}^\infty (-1)^{j-1} e^{-2 j^2 \lambda^2}$. For common significance levels $\alpha$, the asymptotic critical values $\lambda_\alpha$ of $\sqrt{n} D_n$ are $\lambda_{0.10} = 1.224$, $\lambda_{0.05} = 1.358$, and $\lambda_{0.01} = 1.628$ (so $d_{n, 0.05} \approx 1.358 / \sqrt{n}$).

Nikolai Smirnov also extended this framework to the Two-Sample Kolmogorov-Smirnov Problem: given two independent random samples $X_1, \dots, X_{n_1} \stackrel{\text{i.i.d.}}{\sim} F(x)$ and $Y_1, \dots, Y_{n_2} \stackrel{\text{i.i.d.}}{\sim} G(x)$ with empirical distribution functions $F_{n_1}(x)$ and $G_{n_2}(x)$, we test $H_0 : F(x) = G(x)$ for all $x \in \mathbb{R}$ via:

Under $H_0$ when $F = G$ is continuous, all $\binom{n_1 + n_2}{n_1}$ orderings of the combined sample of size $N = n_1 + n_2$ are equally likely, rendering $D_{n_1, n_2}$ strictly distribution-free in finite samples, and as $n_1, n_2 \to \infty$, the scaled statistic converges to the exact same Kolmogorov distribution:

Although the Kolmogorov-Smirnov statistic $D_n = \sup_x |F_n(x) - F_0(x)|$ is intuitively appealing, it suffers from two well-known limitations: (1) it depends only on the single point of maximum vertical deviation rather than integrating deviations across the entire support, and (2) because the variance of the empirical process $\text{Var}(F_n(x)) = F_0(x)[1 - F_0(x)] / n$ reaches its maximum $1/(4n)$ at the median ($F_0(x) = 0.5$) and vanishes at both tails ($F_0(x) \to 0$ or $1$), $D_n$ is most sensitive to shifts near the center of the distribution and exhibits relatively low power against heavy-tailed deviations! To overcome this tail insensitivity, Harald Cramér (1928), Richard von Mises (1931), Nikolai Smirnov (1936), and Theodore W. Anderson and Donald A. Darling (1952) introduced the **Quadratic Empirical Distribution Function Family** based on weighted $L^2$ norms of $F_n(x) - F_0(x)$: $$Q_n = n \int_{-\infty}^\infty \left[ F_n(x) - F_0(x) \right]^2 \psi\left(F_0(x)\right) dF_0(x) \xrightarrow{d} \int_0^1 \mathbb{B}^2(t) \psi(t) \, dt$$ where $\psi(t) \ge 0$ is a non-negative weight function on $(0, 1)$: 1. **Cramér-von Mises Statistic ($\omega^2_n$, with uniform weight $\psi(t) = 1$)**: $$\omega^2_n = n \int_{-\infty}^\infty \left[ F_n(x) - F_0(x) \right]^2 dF_0(x) = \frac{1}{12n} + \sum_{i=1}^n \left[ F_0\left(X_{(i)}\right) - \frac{2i - 1}{2n} \right]^2$$ Notice how simply $\omega^2_n$ is computed from the ordered uniform transforms $U_{(i)} = F_0(X_{(i)})$ by comparing them against the mid-ranks $(2i-1)/(2n)$! By the Karhunen-Loève orthogonal eigenfunction expansion of the Brownian bridge $\mathbb{B}(t) = \sum_{j=1}^\infty \frac{\sqrt{2}}{\pi j} \sin(\pi j t) Z_j$ (where $Z_j \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, 1)$), Parseval's identity yields $\omega^2_n \xrightarrow{d} \int_0^1 \mathbb{B}^2(t) \, dt = \sum_{j=1}^\infty \frac{Z_j^2}{\pi^2 j^2}$. 2. **Anderson-Darling Statistic ($A^2_n$, with inverse-variance tail weight $\psi(t) = \frac{1}{t(1 - t)}$)**: By choosing $\psi(t) = [t(1 - t)]^{-1}$ equal to the reciprocal of the Brownian bridge variance, the Anderson-Darling statistic standardizes the empirical process across all quantiles and places heavy emphasis on the upper and lower tails: $$A^2_n = n \int_{-\infty}^\infty \frac{\left[ F_n(x) - F_0(x) \right]^2}{F_0(x)\left[ 1 - F_0(x) \right]} \, dF_0(x) = -n - \frac{1}{n} \sum_{i=1}^n (2i - 1) \left[ \ln F_0\left(X_{(i)}\right) + \ln\left( 1 - F_0\left(X_{(n + 1 - i)}\right) \right) \right]$$ Under $H_0$, $A^2_n \xrightarrow{d} \int_0^1 \frac{\mathbb{B}^2(t)}{t(1 - t)} \, dt = \sum_{j=1}^\infty \frac{Z_j^2}{j(j + 1)}$. What happens when we test the **composite hypothesis of normality** $H_0 : X_1, \dots, X_n \sim \mathcal{N}(\mu, \sigma^2)$ with unknown mean $\mu$ and unknown variance $\sigma^2$? If a practitioner plugs the sample mean $\bar{X}$ and sample standard deviation $S$ into the K-S statistic: $$\hat{D}_n = \sup_{x \in \mathbb{R}} \left| F_n(x) - \Phi\left( \frac{x - \bar{X}}{S} \right) \right|$$ and compares $\hat{D}_n$ against standard Kolmogorov critical values ($1.358 / \sqrt{n}$ at $\alpha = 0.05$), the test is **severely conservative** (true Type I error drops from $0.05$ to $\approx 0.005$, accompanied by massive loss of power) because $\bar{X}$ and $S$ are fitted to the sample itself, pulling $\Phi((x - \bar{X})/S)$ much closer to $F_n(x)$ than the true population $\Phi((x - \mu)/\sigma)$! - **Lilliefors Test for Normality (1967)**: Hubert Lilliefors recognized that while the distribution of $\hat{D}_n$ is much smaller than $D_n$, by location-scale invariance of $\frac{X_{(i)} - \bar{X}}{S}$, the distribution of $\hat{D}_n$ under $H_0$ is still **completely independent of the unknown $(\mu, \sigma^2)$**! Lilliefors derived the exact adjusted critical values for $\hat{D}_n$ via Monte Carlo simulation: at $\alpha = 0.05$, the Lilliefors critical value is $d_{\text{Lilliefors}, 0.05} \approx \frac{0.886}{\sqrt{n}}$ (more than $35\%$ smaller than Kolmogorov's $\frac{1.358}{\sqrt{n}}$!), and at $\alpha = 0.01$, $d_{\text{Lilliefors}, 0.01} \approx \frac{1.031}{\sqrt{n}}$. (Lilliefors 1969 derived the analogous adjustment for the Exponential distribution with estimated mean $\bar{X}$, where $d_{\text{Exp}, 0.05} \approx \frac{1.06}{\sqrt{n}}$). - **Shapiro-Wilk $W$ Test for Normality (1965)**: Samuel Shapiro and Martin Wilk constructed an analysis-of-variance regression test comparing the square of the Best Linear Unbiased Estimator (BLUE) of $\sigma$ from order statistics against the sample sum of squares: $$W = \frac{\left( \sum_{i=1}^n a_i X_{(i)} \right)^2}{\sum_{i=1}^n (X_i - \bar{X})^2}, \quad \text{where } \mathbf{a} = (a_1, \dots, a_n)^\top = \frac{\mathbf{m}^\top \mathbf{V}^{-1}}{\left( \mathbf{m}^\top \mathbf{V}^{-2} \mathbf{m} \right)^{1/2}}$$ Here $\mathbf{m} = (m_1, \dots, m_n)^\top$ is the vector of expected values of standard normal order statistics $m_i = \mathbb{E}[Z_{(i)}]$ and $\mathbf{V} = \text{Cov}(\mathbf{Z}_{(\cdot)})$ is their $n \times n$ covariance matrix. By the Cauchy-Schwarz inequality, $0 < W \le 1$, with values of $W$ close to $1$ supporting normality and small left-tail values ($W \le w_\alpha$) rejecting normality. Extensive simulation studies (Razali & Wah, 2011; Stephens, 1974) demonstrate that the **Shapiro-Wilk $W$ test and the Anderson-Darling $A^2$ test are the two most powerful omnibus tests of normality** across skewed, leptokurtic, and platykurtic alternatives.

A fundamental distinction between parametric tests (such as the $t$-test or Neyman-Pearson test) and nonparametric EDF tests ($D_n, \omega^2_n, A^2_n$) lies in their global consistency against arbitrary fixed alternatives $F_1 \ne F_0$.

Theorem 6.5 (Universal Consistency of Kolmogorov-Smirnov and Quadratic EDF Tests)

Let $F_0$ be a continuous null CDF, and let $F_1$ be any alternative CDF such that $F_1(x_0) \ne F_0(x_0)$ for at least one point $x_0 \in \mathbb{R}$. Then for any fixed significance level $\alpha \in (0, 1)$, the power of the level-$\alpha$ Kolmogorov-Smirnov test approaches $1$ as $n \to \infty$:

Proof:

Let $\Delta = \sup_{x \in \mathbb{R}} |F_1(x) - F_0(x)| \ge |F_1(x_0) - F_0(x_0)| > 0$. By the reverse triangle inequality in the supremum norm $\|\cdot\|_\infty$:

Multiplying by $\sqrt{n}$ gives $\sqrt{n} D_n \ge \sqrt{n} \Delta - \sqrt{n} \|F_n - F_1\|_\infty$. Under $F_1$, by the Dvoretzky-Kiefer-Wolfowitz inequality, $\sqrt{n} \|F_n - F_1\|_\infty = O_p(1)$ is tight in probability, whereas $\sqrt{n} \Delta \to +\infty$ as $n \to \infty$. Consequently, $\mathbb{P}_{F_1}(\sqrt{n} D_n > \lambda_\alpha) \to 1$.

alongside formal hypothesis tests, visual probability plots—specifically the Quantile-Quantile (Q-Q) Plot graphing the sample order statistics $X_{(i)}$ against the theoretical standard quantiles $F_0^{-1}\left(\frac{i - 0.5}{n}\right)$—provide an indispensable diagnostic companion. Notice that the square root of the Shapiro-Wilk statistic $\sqrt{W}$ is essentially the sample Pearson correlation coefficient between the ordered observations $X_{(i)}$ and the normal scores $m_i$ on the normal Q-Q plot!

Smirnov's Exact Finite-Sample Combinatorial Recursion and Steck's Determinant

In finite samples of size $n$, how are the exact critical values $d_{n, \alpha}$ of the Kolmogorov-Smirnov statistics $D_n^+$ and $D_n$ computed without relying on the large-sample Brownian bridge limit? For the one-sided statistic $D_n^+ = \max_{1 \le i \le n} \left( \frac{i}{n} - U_{(i)} \right)$, Z. W. Birnbaum and Fred H. Tingey (1951) proved an exact closed-form combinatorial formula for the survival probability $\mathbb{P}_{H_0}(D_n^+ > c)$ for any threshold $c \in (0, 1)$ by integrating over the ordered uniform polytope $\left\{ (u_1, \dots, u_n) : 0 < u_1 < \dots < u_n < 1 \text{ and } u_i > \frac{i}{n} - c \right\}$ using Abel's generalization of the binomial theorem:

For the two-sample Kolmogorov-Smirnov test with sample sizes $n_1$ and $n_2$ ($N = n_1 + n_2$), when $n_1 = n_2 = m$, the combined sample path can be represented as a lattice walk from $(0, 0)$ to $(m, m)$ taking $m$ unit horizontal steps ($+1, 0$) corresponding to $X$-observations and $m$ unit vertical steps ($0, +1$) corresponding to $Y$-observations. Under the null hypothesis $H_0 : F = G$, all $\binom{2m}{m}$ lattice paths from $(0, 0)$ to $(m, m)$ are strictly equiprobable! The event $\left\{ D_{m, m} < \frac{c}{m} \right\}$ is equivalent to the geometric condition that the lattice path never touches or crosses the two diagonal boundary lines $y = x + c$ and $y = x - c$. By applying the André Reflection Principle repeatedly across the two diagonal boundaries $y - x = \pm c$, Nikolai Smirnov and B. V. Gnedenko derived the exact finite-sample combinatorial distribution:

This exact reflection series explains both why the finite-sample distribution of $D_{m, m}$ is purely combinatorial and how the alternating series $K(\lambda) = \sum_{j = -\infty}^\infty (-1)^j e^{-2 j^2 \lambda^2}$ emerges naturally in the Stirling asymptotic limit as $m \to \infty$ with $c = \lambda \sqrt{2m}$.

Neyman's Smooth Test of Goodness-of-Fit and Orthogonal Legendre Polynomials

To bridge the gap between omnibus EDF tests (which spread their power thinly across infinite-dimensional space) and parametric likelihood ratio tests (which concentrate power along a single parametric direction), Jerzy Neyman (1937) introduced the Smooth Test of Goodness-of-Fit. Transforming the raw observations via $U_i = F_0(X_i) \in [0, 1]$, testing $H_0 : F = F_0$ is equivalent to testing whether $U_1, \dots, U_n \stackrel{\text{i.i.d.}}{\sim} \text{Uniform}(0, 1)$. Neyman embedded the uniform density $g_0(u) = 1$ on $[0, 1]$ into a smooth $d$-dimensional exponential family generated by the first $d$ orthonormal Legendre polynomials $\pi_1(u), \pi_2(u), \dots, \pi_d(u)$ on $[0, 1]$, where $\int_0^1 \pi_r(u) \, du = 0$ and $\int_0^1 \pi_r(u) \pi_s(u) \, du = \delta_{rs}$:

• $\pi_1(u) = \sqrt{12}\left(u - \frac{1}{2}\right)$ (captures linear shift in location/mean)
• $\pi_2(u) = \sqrt{5}\left[6\left(u - \frac{1}{2}\right)^2 - \frac{1}{2}\right]$ (captures quadratic dispersion/scale change)
• $\pi_3(u) = \sqrt{7}\left[20\left(u - \frac{1}{2}\right)^3 - 3\left(u - \frac{1}{2}\right)\right]$ (captures cubic skewness asymmetry)
• $\pi_4(u) = 3\left[70\left(u - \frac{1}{2}\right)^4 - 15\left(u - \frac{1}{2}\right)^2 + \frac{3}{8}\right]$ (captures quartic kurtosis/tail heaviness)

Under the smooth exponential family density $g(u; \boldsymbol{\theta}) = \exp\left( \sum_{r=1}^d \theta_r \pi_r(u) - \psi(\boldsymbol{\theta}) \right)$, the Rao Score Test statistic for testing $H_0 : \theta_1 = \theta_2 = \dots = \theta_d = 0$ decomposes into the sum of squares of $d$ orthogonal sample components:

By choosing a small order $d \in \{2, 3, 4\}$, Neyman's smooth test achieves substantially higher local Pitman power against low-frequency location, scale, skewness, and kurtosis alternatives than either Pearson's $\chi^2$ test with many bins or the Kolmogorov-Smirnov supremum test, while decomposing the departure from $H_0$ into interpretable orthogonal moment components!

Asymptotic Pitman and Bahadur Efficiency of Goodness-of-Fit Procedures

When selecting among competing goodness-of-fit procedures in applied scientific work, how do Pearson's $\chi^2$ test, the Kolmogorov-Smirnov supremum test $D_n$, the Cramér-von Mises quadratic test $\omega^2_n$, and the Anderson-Darling tail-weighted test $A^2_n$ compare in terms of local Pitman power and nonlocal Bahadur efficiency? Let us examine the fundamental structural trade-offs revealed by asymptotic empirical process theory:

1. Discretization Loss in Pearson's $\chi^2$ Test: When the underlying random variable $X$ is genuinely continuous, partitioning the real line into $k$ intervals discards all intra-interval positional information. Moreover, as the number of cells $k$ increases, the degrees of freedom $k - 1$ of the limiting central $\chi^2_{k-1}$ distribution grow, which increases the null variance $2(k - 1)$ of the test statistic and dilutes the noncentrality parameter $\lambda = n \sum_{j=1}^k \frac{(p_{j1} - p_{j0})^2}{p_{j0}}$ across more degrees of freedom. Consequently, against smooth location or scale shifts of a continuous distribution, Pearson's $\chi^2$ test has strictly zero Pitman Asymptotic Relative Efficiency ($\text{ARE} = 0$) relative to parametric score tests unless the bin boundaries are specifically optimized, and systematically trails Empirical Distribution Function tests such as $\omega^2_n$ and $A^2_n$.

2. Supremum ($L^\infty$) versus Integrated Quadratic ($L^2$) EDF Statistics: Because the Kolmogorov-Smirnov statistic $D_n = \sup_{x \in \mathbb{R}} |F_n(x) - F_0(x)|$ evaluates only the single largest pointwise vertical discrepancy, its sensitivity is governed by the supremum of the deterministic shift function $\delta(x) = F_1(x) - F_0(x)$. In contrast, the Cramér-von Mises statistic $\omega^2_n$ integrates $[\delta(x)]^2$ over the entire support with respect to $dF_0(x)$, accumulating smooth global departures across every quantile. By the Karhunen-Loève spectral decomposition of the Brownian bridge $\mathbb{B}(t) = \sum_{j=1}^\infty \frac{\sqrt{2}}{\pi j} \sin(\pi j t) Z_j$, the first principal component $Z_1$ (which corresponds to a global location shift $\sin(\pi t)$) receives weight $\frac{1}{\pi^2} \approx 0.1013$ out of the total variance $\sum_{j=1}^\infty \frac{1}{\pi^2 j^4} = \frac{1}{90}$ in $\omega^2_n$, accounting for over $91\%$ of the total null variance of the Cramér-von Mises statistic! This spectral concentration explains why $\omega^2_n$ and $A^2_n$ achieve substantially higher statistical power than $D_n$ against smooth shifts in mean, variance, skewness, and tail thickness.

---