Stein's method

Stein's method is a way to bound the distance between two probability distributions using a specific probability metric.

Let the metric be given in the form

${\displaystyle (1.1)\quad d(P,Q)=\sup _{h\in {\mathcal {H}}}\left|\int hdP-\int hdQ\right|=\sup _{h\in {\mathcal {H}}}\left|Eh(W)-Eh(Y)\right|}$

Here, ${\displaystyle P}$ and ${\displaystyle Q}$ are probability measures on a measurable space ${\displaystyle {\mathcal {X}}}$, ${\displaystyle W}$ and ${\displaystyle Y}$ are random variables with distribution ${\displaystyle P}$ and ${\displaystyle Q}$ respectively, ${\displaystyle E}$ is the usual expectation operator and ${\displaystyle {\mathcal {H}}}$ is a set of functions from ${\displaystyle {\mathcal {X}}}$ to the set of real numbers. Set ${\displaystyle {\mathcal {H}}}$ has to be large enough, so that the above definition indeed yields a metric.

Important examples are the total variation metric, where we let ${\displaystyle {\mathcal {H}}}$ consist of all the indicator functions of measurable sets, the Kolmogorov (uniform) metric for probability measures on the real numbers, where we consider all the half-line indicator functions, and the Lipschitz (first order Wasserstein; Kantorovich) metric, where the underlying space is itself a metric space and we take the set ${\displaystyle {\mathcal {H}}}$ to be all Lipschitz-continuous functions with Lipschitz-constant 1. However, note that not every metric can be represented in the form (1.1).

In what follows ${\displaystyle P}$ is a complicated distribution (e.g., the distribution of a sum of dependent random variables), which we want to approximate by a much simpler and tractable distribution ${\displaystyle Q}$ (e.g., the standard normal distribution).

We assume now that the distribution ${\displaystyle Q}$ is a fixed distribution; in what follows we shall in particular consider the case where ${\displaystyle Q}$ is the standard normal distribution, which serves as a classical example.

First of all, we need an operator ${\displaystyle {\mathcal {A}}}$, which acts on functions ${\displaystyle f}$ from ${\displaystyle {\mathcal {X}}}$ to the set of real numbers and 'characterizes' distribution ${\displaystyle Q}$ in the sense that the following equivalence holds:

${\displaystyle (2.1)\quad E({\mathcal {A}}f)(Y)=0{\text{ for all }}f\quad \iff \quad Y{\text{ has distribution }}Q.}$

We call such an operator the Stein operator.

For the standard normal distribution, Stein's lemma yields such an operator:

${\displaystyle (2.2)\quad E\left(f'(Y)-Yf(Y)\right)=0{\text{ for all }}f\in C_{b}^{1}\quad \iff \quad Y{\text{ has standard normal distribution.}}}$

Thus, we can take

${\displaystyle (2.3)\quad ({\mathcal {A}}f)(x)=f'(x)-xf(x).}$

There are in general infinitely many such operators and it still remains an open question, which one to choose. However, it seems that for many distributions there is a particular good one, like (2.3) for the normal distribution.

There are different ways to find Stein operators.[4]

${\displaystyle P}$ is close to ${\displaystyle Q}$ with respect to ${\displaystyle d}$ if the difference of expectations in (1.1) is close to 0. We hope now that the operator ${\displaystyle {\mathcal {A}}}$ exhibits the same behavior: if ${\displaystyle P=Q}$ then ${\displaystyle E({\mathcal {A}}f)(W)=0}$, and hopefully if ${\displaystyle P\approx Q}$ we have ${\displaystyle E({\mathcal {A}}f)(W)\approx 0}$.

It is usually possible to define a function ${\displaystyle f=f_{h}}$ such that

${\displaystyle (3.1)\quad ({\mathcal {A}}f)(x)=h(x)-E[h(Y)]\qquad {\text{ for all }}x.}$

We call (3.1) the Stein equation. Replacing ${\displaystyle x}$ by ${\displaystyle W}$ and taking expectation with respect to ${\displaystyle W}$, we get

${\displaystyle (3.2)\quad E({\mathcal {A}}f)(W)=E[h(W)]-E[h(Y)].}$

Now all the effort is worthwhile only if the left-hand side of (3.2) is easier to bound than the right hand side. This is, surprisingly, often the case.

If ${\displaystyle Q}$ is the standard normal distribution and we use (2.3), then the corresponding Stein equation is

${\displaystyle (3.3)\quad f'(x)-xf(x)=h(x)-E[h(Y)]\qquad {\text{for all }}x.}$

If probability distribution Q has an absolutely continuous (with respect to the Lebesgue measure) density q, then[4]

${\displaystyle (3.4)\quad ({\mathcal {A}}f)(x)=f'(x)+f(x)q'(x)/q(x).}$

Analytic methods. Equation (3.3) can be easily solved explicitly:

${\displaystyle (4.1)\quad f(x)=e^{x^{2}/2}\int _{-\infty }^{x}[h(s)-Eh(Y)]e^{-s^{2}/2}ds.}$

Generator method. If ${\displaystyle {\mathcal {A}}}$ is the generator of a Markov process ${\displaystyle (Z_{t})_{t\geq 0}}$ (see Barbour (1988), Götze (1991)), then the solution to (3.2) is

${\displaystyle (4.2)\quad f(x)=-\int _{0}^{\infty }[E^{x}h(Z_{t})-Eh(Y)]dt,}$

where ${\displaystyle E^{x}}$ denotes expectation with respect to the process ${\displaystyle Z}$ being started in ${\displaystyle x}$. However, one still has to prove that the solution (4.2) exists for all desired functions ${\displaystyle h\in {\mathcal {H}}}$.

Usually, one tries to give bounds on ${\displaystyle f}$ and its derivatives (or differences) in terms of ${\displaystyle h}$ and its derivatives (or differences), that is, inequalities of the form

${\displaystyle (5.1)\quad \|D^{k}f\|\leq C_{k,l}\|D^{l}h\|,}$

for some specific ${\displaystyle k,l=0,1,2,\dots }$ (typically ${\displaystyle k\geq l}$ or ${\displaystyle k\geq l-1}$, respectively, depending on the form of the Stein operator), where often ${\displaystyle \|\cdot \|}$ is the supremum norm. Here, ${\displaystyle D^{k}}$ denotes the differential operator, but in discrete settings it usually refers to a difference operator. The constants ${\displaystyle C_{k,l}}$ may contain the parameters of the distribution ${\displaystyle Q}$. If there are any, they are often referred to as Stein factors.

In the case of (4.1) one can prove for the supremum norm that

${\displaystyle (5.2)\quad \|f\|_{\infty }\leq \min\{{\sqrt {\pi /2}}\|h\|_{\infty },2\|h'\|_{\infty }\},\quad \|f'\|_{\infty }\leq \min\{2\|h\|_{\infty },4\|h'\|_{\infty }\},\quad \|f''\|_{\infty }\leq 2\|h'\|_{\infty },}$

where the last bound is of course only applicable if ${\displaystyle h}$ is differentiable (or at least Lipschitz-continuous, which, for example, is not the case if we regard the total variation metric or the Kolmogorov metric!). As the standard normal distribution has no extra parameters, in this specific case the constants are free of additional parameters.

If we have bounds in the general form (5.1), we usually are able to treat many probability metrics together. One can often start with the next step below, if bounds of the form (5.1) are already available (which is the case for many distributions).

We are now in a position to bound the left hand side of (3.1). As this step heavily depends on the form of the Stein operator, we directly regard the case of the standard normal distribution.

At this point we could directly plug in random variable ${\displaystyle W}$, which we want to approximate, and try to find upper bounds. However, it is often fruitful to formulate a more general theorem. Consider here the case of local dependence.

Assume that ${\displaystyle W=\sum _{i=1}^{n}X_{i}}$ is a sum of random variables such that the ${\displaystyle E[W]=0}$ and variance ${\displaystyle \operatorname {var} [W]=1}$. Assume that, for every ${\displaystyle i=1,\dots ,n}$, there is a set ${\displaystyle A_{i}\subset \{1,2,\dots ,n\}}$, such that ${\displaystyle X_{i}}$ is independent of all the random variables ${\displaystyle X_{j}}$ with ${\displaystyle j\not \in A_{i}}$. We call this set the 'neighborhood' of ${\displaystyle X_{i}}$. Likewise let ${\displaystyle B_{i}\subset \{1,2,\dots ,n\}}$ be a set such that all ${\displaystyle X_{j}}$ with ${\displaystyle j\in A_{i}}$ are independent of all ${\displaystyle X_{k}}$, ${\displaystyle k\not \in B_{i}}$. We can think of ${\displaystyle B_{i}}$ as the neighbors in the neighborhood of ${\displaystyle X_{i}}$, a second-order neighborhood, so to speak. For a set ${\displaystyle A\subset \{1,2,\dots ,n\}}$ define now the sum ${\displaystyle X_{A}:=\sum _{j\in A}X_{j}}$.

Using Taylor expansion, it is possible to prove that

${\displaystyle (6.1)\quad \left|E(f'(W)-Wf(W))\right|\leq \|f''\|_{\infty }\sum _{i=1}^{n}\left({\frac {1}{2}}E|X_{i}X_{A_{i}}^{2}|+E|X_{i}X_{A_{i}}X_{B_{i}\setminus A_{i}}|+E|X_{i}X_{A_{i}}|E|X_{B_{i}}|\right)}$

Note that, if we follow this line of argument, we can bound (1.1) only for functions where ${\displaystyle \|h'\|_{\infty }}$ is bounded because of the third inequality of (5.2) (and in fact, if ${\displaystyle h}$ has discontinuities, so will ${\displaystyle f''}$). To obtain a bound similar to (6.1) which contains only the expressions ${\displaystyle \|f\|_{\infty }}$ and ${\displaystyle \|f'\|_{\infty }}$, the argument is much more involved and the result is not as simple as (6.1); however, it can be done.

Theorem A. If ${\displaystyle W}$ is as described above, we have for the Lipschitz metric ${\displaystyle d_{W}}$ that

${\displaystyle (6.2)\quad d_{W}({\mathcal {L}}(W),N(0,1))\leq 2\sum _{i=1}^{n}\left({\frac {1}{2}}E|X_{i}X_{A_{i}}^{2}|+E|X_{i}X_{A_{i}}X_{B_{i}\setminus A_{i}}|+E|X_{i}X_{A_{i}}|E|X_{B_{i}}|\right).}$

Proof. Recall that the Lipschitz metric is of the form (1.1) where the functions ${\displaystyle h}$ are Lipschitz-continuous with Lipschitz-constant 1, thus ${\displaystyle \|h'\|\leq 1}$. Combining this with (6.1) and the last bound in (5.2) proves the theorem.

Thus, roughly speaking, we have proved that, to calculate the Lipschitz-distance between a ${\displaystyle W}$ with local dependence structure and a standard normal distribution, we only need to know the third moments of ${\displaystyle X_{i}}$ and the size of the neighborhoods ${\displaystyle A_{i}}$ and ${\displaystyle B_{i}}$.

We can treat the case of sums of independent and identically distributed random variables with Theorem A.

Assume that ${\displaystyle EX_{i}=0}$, ${\displaystyle varX_{i}=1}$ and ${\displaystyle W=n^{-1/2}\sum X_{i}}$. We can take ${\displaystyle A_{i}=B_{i}=\{i\}}$. From Theorem A we obtain that

${\displaystyle (7.1)\quad d_{W}({\mathcal {L}}(W),N(0,1))\leq {\frac {5E|X_{1}|^{3}}{n^{1/2}}}.}$

For sums of random variables another approach related to Steins Method is known as the zero bias transform.