Epps–Pulley Statistic: Normality Test

• The Epps–Pulley statistic is a test of normality based on the weighted L²-distance between the empirical characteristic function of standardized data and that of a standard normal distribution.
• It naturally extends to multivariate settings via affine-invariant standardizations and is equivalent to a biased Gaussian-kernel maximum mean discrepancy under simple transformations.
• Its well-understood asymptotic behavior, tunable smoothing parameter β, and spectral representation enable efficient p-value computation and diagnostic use in diverse modeling applications.

The Epps–Pulley statistic is a test of normality constructed as a weighted $L^2$-distance between the empirical characteristic function of standardized data and that of the standard normal distribution, with a Gaussian weight parameterized by a bandwidth parameter $\beta>0$. It extends naturally to multivariate (affine-invariant, universally consistent) settings and is closely related, under simple transformations, to the biased Gaussian-kernel maximum mean discrepancy (MMD) between an empirical distribution and the standard normal. The Epps–Pulley statistic is distribution-free under normality and has well-understood asymptotic behavior, efficiency, and practical performance across diverse alternatives to normality (Ebner et al., 2021, Ebner et al., 2021, Rustamov, 2019).

1. Definition and Forms of the Epps–Pulley Statistic

For a univariate i.i.d. sample $X_1,...,X_n$ with unknown mean $\mu$ and variance $\sigma^2$, let $\overline X_n = n^{-1}\sum_{j=1}^n X_j$ and $$S_n^2 = n^{-1} \sum_{j=1}^n (X_j - \overline X_n)^2$$. The scaled residuals are $Y_{n,j} = (X_j - \overline X_n)/S_n$. The empirical characteristic function is $\psi_n(t) = n^{-1} \sum_{j=1}^n e^{i t Y_{n,j}}$, and the Gaussian weight is $$\varphi_\beta(t) = (2\pi \beta^2)^{-1/2} \exp(-t^2/(2\beta^2))$$.

The univariate Epps–Pulley statistic is

$$T_{n,\beta} = n \int_{-\infty}^\infty |\psi_n(t) - e^{-t^2/2}|^2 \varphi_\beta(t) \, dt$$

which can be explicitly computed as

$$T_{n,\beta} = \frac{1}{n}\sum_{j,k=1}^n \exp\left( -\frac{\beta^2}{2}(Y_{n,j} - Y_{n,k})^2 \right) - \frac{2}{\sqrt{1+\beta^2}} \sum_{j=1}^n \exp\left( -\frac{\beta^2 Y_{n,j}^2}{2(1+\beta^2)} \right) + \frac{n}{\sqrt{1+2\beta^2}}$$

For the multivariate case, if $X_1, ..., X_n \in \mathbb{R}^d$ have sample mean $\overline X_n$ and covariance matrix $\Sigma_n$, then standardized residuals $Y_{n,j} = \Sigma_n^{-1/2}(X_j - \overline X_n)$ are used and the integration occurs over $\mathbb{R}^d$ with $$\varphi_{\beta,d}(t) = (2\pi \beta^2)^{-d/2} \exp( -\|t\|^2/2\beta^2 )$$(Ebner et al., 2021, Ebner et al., 2021). This generalization preserves affine invariance and universal consistency.

2. Connection to Maximum Mean Discrepancy (MMD)

The Epps–Pulley statistic coincides, up to notation and centering, with the biased Gaussian-kernel MMD between an empirical distribution $Q_n$ and the standard normal $\mathcal N_d(0,I)$:

$$\mathrm{MMD}_b^2(Q_n; \gamma) = \left(\frac{\gamma^2}{2+\gamma^2}\right)^{d/2} - \frac{2}{n} \left(\frac{\gamma^2}{1+\gamma^2}\right)^{d/2} \sum_{i=1}^n \exp\left( -\frac{\|z_i\|^2}{2(1+\gamma^2)} \right) + \frac{1}{n^2} \sum_{i,j} \exp\left( -\frac{\|z_i-z_j\|^2}{2\gamma^2} \right)$$

with the correspondence $\gamma = 1/\beta$ (Rustamov, 2019). The Epps–Pulley (or equivalently, the BHEP) statistic thus inherits the omnibus consistency and interpretability of the MMD, where the smoothing parameter $\beta$ (equivalently, kernel width $\gamma$) modulates sensitivity to different departures from normality.

3. Asymptotic Null Distribution and Eigenstructure

Under $H_0\!:\, X_j \sim \mathcal N(\mu, \sigma^2)$, the Epps–Pulley statistic $T_{n,\beta}$ is asymptotically distribution-free. As $n \rightarrow \infty$,

$$T_{n,\beta} \xrightarrow{d} T_\beta := \int_{\mathbb R} Z(t)^2 \varphi_\beta(t) dt$$

where $Z(\cdot)$ is a centered Gaussian process with covariance kernel

$$K(s, t) = \exp\left( -\frac{(s-t)^2}{2} \right) - [1 + st + (st)^2/2] \exp\left( -\frac{s^2 + t^2}{2} \right)$$

(Ebner et al., 2021, Ebner et al., 2021). The limiting law admits the series representation

$$T_\beta \stackrel{d}{=} \sum_{j=1}^\infty \lambda_j(\beta) N_j^2$$

where $N_j$ are i.i.d. $N(0,1)$ and $\{\lambda_j(\beta)\}$ are the positive eigenvalues of the integral operator on $L^2(\mathbb R, \varphi_\beta)$ defined by $K$. Numerical values for $\beta=1$ include $$\lambda_1 \approx 0.0743, \lambda_2 \approx 0.0448, ...$$(Ebner et al., 2021). This spectral representation underpins tail approximations, p-value computations, and Bahadur slopes.

4. Bahadur Efficiency and Comparative Power

The Bahadur efficiency of the Epps–Pulley statistic against local alternatives $G(\cdot;\theta)$ indexed by $\theta$ is characterized by two quantities:

• $b_{T_\beta}(\theta)$, the asymptotic mean under $G(\cdot; \theta)$,
• $\lambda_1(\beta)$, the leading eigenvalue from the null distribution.

Specifically, for small $\theta$,

$$b_{T_\beta}(\theta) \approx \Delta(\beta) \theta^2, \qquad c_{T_\beta}^*(\theta) = \frac{b_{T_\beta}(\theta)^2}{\lambda_1(\beta)}, \qquad \ell(\beta) = \frac{\Delta(\beta)^2}{\lambda_1(\beta)}$$

(Ebner et al., 2021). Table 1 in (Ebner et al., 2021) demonstrates that for $\beta$ in the range $0.5 - 1$, the efficiency of the Epps–Pulley test is competitive with or generally superior to EDF-based tests (e.g., Shapiro–Wilk, Anderson–Darling) for certain alternatives, especially those exhibiting skewness (Ley–Paindaveine alternatives):

β	Lehmann	1st Ley–P.	2nd Ley–P.	Contam. N(1,1)	Contam. N(.5,1)	Contam. N(0,.5)
0.25	.996	.947	.824	.760	.945	.084
0.5	.895	.944	.872	.649	.824	.267
0.75	.854	.998	.986	.592	.766	.474
1	.743	.937	.981	.499	.654	.587

Values are efficiencies relative to the (infeasible) likelihood-ratio test (Ebner et al., 2021). This suggests optimal power for $\beta$ near unity across a range of alternatives.

5. Practical Guidelines for Tuning and Implementation

The smoothing parameter $\beta$ modulates the Epps–Pulley test’s sensitivity:

• Small $\beta$ ($0.25$–$0.5$): Emphasizes high-frequency deviations in the characteristic function, suitable for location-type alternatives and mild skew.
• Intermediate $\beta$ ($0.75$–$1$): Effective across both skew and scale-contamination alternatives; $\beta=1$ aligns with classical BHEP settings.
• Large $\beta$ ($2$–$5$): Favours low-frequency discrepancies, enhancing sensitivity to scale or heavy-tailed alternatives at the cost of skew-detection (Ebner et al., 2021, Rustamov, 2019).

Implementation requires selection or cross-validation of $\beta$, or adoption of empirically motivated heuristics (e.g., Henze–Zirkler’s $\gamma_{d,n}$ for high-dimensional settings). Under the composite null, centering and whitening are performed prior to evaluation (Rustamov, 2019).

Monte Carlo can be used to estimate finite-sample thresholds or compute the Bahadur efficiency via plug-in estimates of $b_{T_\beta}(\theta)$ and $\lambda_1(\beta)$. In practice, standardized versions using explicit null mean/variance formulas allow rapid computation and interpretation (Rustamov, 2019).

6. Comparative and Theoretical Properties

The Epps–Pulley statistic is affine-invariant and, in the multivariate case, universally consistent against alternatives with finite second moments (Ebner et al., 2021, Rustamov, 2019). Compared to Shapiro–Wilk and Anderson–Darling, its key distinctions are:

• Direct applicability in arbitrary dimensions,
• Closed-form critical value calculations via eigen-expansions,
• Explicit power modulated by the smoothing parameter,
• One-to-one equivalence (up to bias) with the widely used MMD, situating it within reproducing kernel Hilbert space (RKHS) theoretical frameworks.

Under the simple null $X_i \sim \mathcal N_d(0, I)$, standardized Epps–Pulley values are approximately $N(0,1)$; under the composite null, centering/whitening preserves distribution-freeness up to small sample effects (Rustamov, 2019).

A plausible implication is that the flexibility of the smoothing parameter $\beta$ allows the Epps–Pulley statistic to target classes of alternatives of interest, further justifying its use in modern model criticism and generative modeling settings.

7. Computational Considerations and Extensions

The computation of the Epps–Pulley statistic is explicit and involves only matrix operations and exponential functions. For limiting distribution evaluation, the eigenvalues $\lambda_j(\beta)$ can be computed numerically by truncating the spectral expansion, with software implementation benefiting from the work of Ebner and Henze (Ebner et al., 2021). For p-values, Imhof’s method or cumulant-matching to the Pearson system is effective.

In the context of generative models such as Wasserstein Auto-Encoders, batch-averaged and standardized forms of the Epps–Pulley (MMD) statistic are employed for both training objectives and diagnostics, with the correspondence between kernel width $\gamma$ and smoothing $\beta$ guiding model selection procedures (Rustamov, 2019).

---

References:

(Ebner et al., 2021) "Bahadur efficiencies of the Epps--Pulley test for normality" (Ebner et al., 2021) "On the eigenvalues associated with the limit null distribution of the Epps-Pulley test of normality" (Rustamov, 2019) "Closed-form Expressions for Maximum Mean Discrepancy with Applications to Wasserstein Auto-Encoders"