Stein’s Method
A normal approximation can be tested through an integration-by-parts identity. Stein’s method estimates the error in that identity and then converts the estimate into a distance between distributions. We begin with independent sums, where removing one summand gives the necessary cancellation, and then extend the argument to local dependence. For the preceding discussion of dependence, see mixing conditions and sample means.
A useful starting point is the following characterization of a random variable \(W\) having the standard normal distribution. A random variable \(W\) satisfies
for all “smooth” functions \(f\) if and only if \(W\) has the standard normal distribution. This is a consequence of Gaussian integration by parts for bounded continuously differentiable \(f\) with bounded derivative. In some sense, \(W\) is a standard normal random variable if \(\mathbb{E}[f'(W)] - \mathbb{E}[W f(W)] \approx 0\) in a suitable sense. This is possible due to the following heuristic: let \(X_1,\dots,X_n\) be i.i.d. random variables with mean zero and variance one. Define
Seen this way, \(X_i \perp \!\!\! \perp W_i\) by construction. We therefore can expect, because \(\mathbb{E}[X_i f(W_i)] = \mathbb{E}[X_i]\mathbb{E}f(W_i) = 0\),
by Taylor’s theorem. Therefore,
where the first approximate equality in the second line follows from Taylor’s theorem, and the second approximate equality in the second line follows from the law of large numbers: \(\frac1n \sum_{1 \leq i \leq n} X_i^2 \approx 1\). In fact, a preliminary formalization of this heuristic is the following result that essentially reverses the Gaussian integration by parts.
Proposition 1. If \(Z\sim \mathcal{N}(0,1)\), then \(\mathbb{E}f^{\prime}(Z)=\mathbb{E}[Zf(Z)]\) for every bounded continuously differentiable \(f\) with bounded derivative. Conversely, if \(\mathbb{E}|W|<\infty\) and this identity holds with \(W\) in place of \(Z\) for all such tests, then \(W\sim \mathcal{N}(0,1)\).
Proof.
For the forward direction integrate against the Gaussian density, whose derivative is \(-x\varphi(x)\); boundedness makes the boundary term zero. For the converse, apply the identity to the real and imaginary parts of \(f(w)=e^{itw}\). The characteristic function \(\chi(t)=Ee^{itW}\) is differentiable because \(\mathbb{E}|W|<\infty\), and satisfies \(\chi^{\prime}(t)=-t\chi(t)\) with \(\chi(0)=1\). Therefore \(\chi(t)=e^{-t^2/2}\), which characterizes the standard normal law.
◻
To turn the heuristic into a bound on distribution functions, choose a differential equation whose right-hand side is the difference between an indicator and its Gaussian expectation.
Lemma 2. For each \(x\in\mathbb R\), the following equation has a unique bounded, locally absolutely continuous solution, with the equation understood almost everywhere in \(w\):
is given by the formula
To verify boundedness and uniqueness in Lemma 2, multiplying by \(e^{-w^2/2}\) gives an absolutely continuous primitive. The integral of the centered right-hand side over the whole real line is zero, so the lower-tail expression can equally be written as the negative of the upper-tail integral. For \(w>0\),
since \(t/w\ge1\) inside the integral and \(\int_w^\infty t e^{-t^2/2}\,dt=e^{-w^2/2}\). The reflected inequality handles the negative tail; continuity handles any fixed compact interval. This proves boundedness. The difference of two solutions satisfies \(u'=wu\) almost everywhere, hence \(u(w)=c e^{w^2/2}\). Boundedness forces \(c=0\). Evaluating the primitive separately on \(w\le x\) and \(w>x\) yields the two displayed formulas.
The solution need not have a classical derivative at \(w=x\). For the next identity, define the value of \(f_x^{\prime}\) at that point by the right-hand side of the Stein equation plus \(xf_x(x)\). This derivative version makes the identity valid even when \(W\) has an atom at \(x\). Taking expectations gives
so the problem of proving (the “finite sample validity” of) the approximation of \(W\) using the standard normal distribution becomes bounding \(|\mathbb{E}[f'(W) - Wf(W)]|\) from above. The choice of test functions determines the distance between distributions. Let \(\mathcal H\) denote some function class, and we are interested in answering under what conditions on \(\mathcal H\) is \(\mathbb{E}h(X) \approx \mathbb{E}h(Y)\) for random variables \(X\) and \(Y\). In particular, consider the metric
For a random variable \(W\) and \(Z \sim \mathcal{N}(0, 1)\), in order to compare \(\mathbb{E}h(W)\) with \(\mathbb{E}h(Z)\), fix a function \(h\) with \(\mathbb{E}|h(Z)| < \infty\). We can adapt the construction in Lemma 2 to obtain a corresponding solution to the differential equation
which is known as the Stein equation. In particular, we have the following result.
Proposition 3. The solution to (1) is written as
Furthermore, for all bounded functions \(h\),
If \(h\) is absolutely continuous with \(\|h^{\prime}\|_\infty<\infty\), boundedness of \(h\) is unnecessary, and we have
Multiplying the Stein equation by \(e^{-w^2/2}\) and integrating verifies the displayed solution; its centering is \(\mathbb{E}h(Z)\). The derivative estimates in Proposition 3, with essential supremum norms and finite \(\|h^{\prime}\|_\infty\) in the Lipschitz case, are Ross (2011), Lemma 2.5. We use these exact estimates as an external lemma. They give the following bound for independent sums.
Theorem 4. Let \(X_1,\ldots,X_n\) be independent, centered, with variance one and finite absolute third moments. Put \(W=n^{-1/2}\sum_iX_i\). If \(f\) has a Lipschitz derivative with constant \(C\), then
Proof.
Set \(W_i=W-X_i/\sqrt n\). Independence and centering give \(\mathbb{E}[X_if(W_i)]=0\). Taylor expansion at \(W_i\) and \(\mathbb{E}X_i^2=1\) imply
The remaining difference from \(\mathbb{E}f^{\prime}(W)\) is at most \(Cn^{-3/2}\sum_i\mathbb{E}|X_i|\). Since \(\mathbb{E}|X_i|\leq1\leq \mathbb{E}|X_i|^3\), adding these bounds proves the claim.
◻
The Taylor estimate used in Theorem 4 follows from an integral, so it remains valid when the second derivative exists only almost everywhere:
Apply this with \(v=W_i\), \(a=X_i/\sqrt n\), multiply by \(X_i/\sqrt n\), and sum. Independence factors \(E[X_i^2 f'(W_i)]=E f'(W_i)\). The remaining replacement of \(f'(W_i)\) by \(f'(W)\) costs at most \(C|X_i|/\sqrt n\) for each term. Finally, Cauchy–Schwarz gives \(E|X_i|\le1\), and monotonicity of moment norms gives \(E|X_i|^3\ge1\). These are precisely the two contributions \(C/2\) and \(C\) to the stated constant.
For distribution functions, one seeks a bound uniform in the threshold. The following Berry–Esseen estimate gives such a bound under a third-moment assumption.
Theorem 5. Assume \(X_1, \dots , X_n\) are i.i.d. with \(\mathbb{E}X_i = 0\), \(\mathop{\mathrm{Var}}(X_i) = 1\), and \(\mathbb{E}|X_i|^3 < \infty\). Let \(W_n = \frac{1}{\sqrt n} \sum_{1 \leq i \leq n} X_i\). Then,
The numerical constant in Theorem 5 is an external Berry–Esseen bound, not a consequence of the smooth estimate in Theorem 4. It follows, for example, from the stronger constant \(0.4748\) in Shevtsova (2011), Corollary 1, under the same i.i.d., centered, unit-variance and finite-third-moment assumptions. Identical distribution is unnecessary for an \(n^{-1/2}\) rate when suitable moments are bounded uniformly. We next work with \(W_1(W,Z)=\sup_{\operatorname{Lip}(h)\leq1}|\mathbb{E}h(W)-\mathbb{E}h(Z)|\), the Wasserstein distance for integrable variables.
Theorem 6. Let \(X_1, \dots, X_n\) be independent random variables with \(\mathbb{E}X_i = 0\), \(\mathbb{E}X_i^2 = 1\), and \(\mathbb{E}X_i^4 < \infty\). Let \(W = \frac{1}{\sqrt n} \sum_{1 \leq i \leq n} X_i\), and \(Z \sim \mathcal{N}(0,1)\). Then,
The Wasserstein bound reduces to a supremum of Stein discrepancies over functions with controlled derivatives:
Theorem 7. Let \(\mathbb{E}|W|<\infty\) and \(Z\sim \mathcal{N}(0,1)\). Let \(\mathcal H\) consist of bounded continuously differentiable functions satisfying \(\|f\|_\infty\leq2\), \(\|f^{\prime}\|_\infty\leq\sqrt{2/\pi}\), and \(\operatorname{Lip}(f^{\prime})\leq2\). Then
For each 1-Lipschitz test \(h\), Proposition 3 supplies a solution in \(\mathcal H\). Substitution in the Stein equation proves the bound; these are the bounds and function class in Ross, Theorem 3.1.
Proof.
Proof of Theorem 6. Independence gives the exact identity
Taylor’s remainder bound, this time with the derivative evaluated at \(W\), gives
The variance term is bounded by
Apply Theorem 7. If the third and fourth moments are bounded uniformly in \(i\), both terms are \(O(n^{-1/2})\), even without identical distributions.
◻
We can also obtain similar results if we use dependency neighbourhoods. For the bound below, assume that \((X_1,\ldots,X_n)\) has an undirected dependency graph: whenever two disjoint sets of vertices have no edge between them, the corresponding collections of random variables are independent. Let \(N_i\) be the closed neighborhood of \(i\), including \(i\) itself. In particular, \(X_i\) is independent of \(\{X_j:j\notin N_i\}\). The graph assumption also gives the independence between pairs used in the variance calculation below.
Theorem 8. Let \(X_1, \ldots, X_n\) be random variables with \(\mathbb{E}X_i^4 < \infty\), \(\mathbb{E}X_i = 0\), \(\sigma^2= \mathop{\mathrm{Var}}\left( \sum_{1 \leq i \leq n} X_i \right)>0\), and define \(W = \frac{1}{\sigma} \sum_{1 \leq i \leq n} X_i\). Let the collection \((X_1, \ldots, X_n)\) have dependency neighbourhoods \(\{N_i\}_{1 \leq i \leq n}\), with \(D \triangleq \max_{1\leq i \leq n} |N_i|\). Then, for \(Z \sim \mathcal{N}(0,1)\),
Proof.
For \(f\in\mathcal H\) put \(W_i=\sigma^{-1}\sum_{j\notin N_i}X_j\) and \(Q=\sum_i\sum_{j\in N_i}X_iX_j\). The graph assumption gives \(\mathbb{E}[X_if(W_i)]=0\) and \(EQ=\sigma^2\). The same expansion as above yields
The arithmetic–geometric mean inequality bounds each cubic product by one third of the sum of its three absolute cubes. Since the graph is undirected and each closed neighborhood has at most \(D\) vertices, each resulting moment occurs at most \(D^2\) times. This proves the first term in the theorem.
The cubic count can be checked slot by slot. A fixed \(i\) has at most \(D^2\) choices of \((j,k)\). For a fixed middle index \(j\), there are at most \(D\) possible \(i\)'s because \(j\in N_i\) is equivalent to \(i\in N_j\); each allows at most \(D\) choices of \(k\). The last slot has the same bound. Therefore
Moreover, \(EQ=\sigma^2\) because in the full expansion of \(E(\sum_iX_i)^2\), every term with \(j\notin N_i\) vanishes by centering and independence. This verifies the centering of the quadratic error before its variance is bounded.
For the variance, index \(Q\) by ordered pairs \(e=(i,j),\ j\in N_i\), writing \(A_e=X_iX_j\). The covariance of \(A_e\) and \(A_f\) vanishes unless the two pairs overlap or have an edge between them. For fixed \(e=(i,j)\), every such \(f=(k,\ell)\) has at least one endpoint in \(N_i\cup N_j\). There are at most \(4D^2\) ordered pairs of this kind. Using \(|\operatorname{Cov}(A_e,A_f)|\leq(\operatorname{Var}(A_e)+\operatorname{Var}(A_f))/2\) and symmetry of this relation,
The last estimate is weaker than the one just proved, but gives exactly the stated coefficient \(\sqrt{28/\pi}\). Taking the supremum in Theorem 7 completes the proof.
◻