Stein's method for the Beta distribution and the Pólya-Eggenberger Urn

Introduction

The classical Pólya-Eggenberger urn at time zero contains $\alpha\geq 1$ white and $\beta\geq 1$ black balls, and at every positive integer time a ball is chosen uniformly from the urn, independently of the past, and replaced along with $m\geq 1$ additional balls of the same color. With ${\cal L}$ indicating distribution, or law, and $\rightarrow_{d}$ indicating convergence in distribution, it is well known, see for instance, that if $S_{n}=S_{n}^{\alpha,\beta,m}$ is the number of white balls drawn from the urn by time $n=0,1,2,\ldots$ then as $n\rightarrow\infty$

Here, for positive real numbers $\alpha$ and $\beta$ we let ${\cal B}(\alpha,\beta)$ denote the Beta distribution having density

where $B(\alpha,\beta)=\Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)$ is the Beta function as expressed in terms of the Gamma function $\Gamma(x)$.

Using Stein’s method we derive an order $O(1/n)$ bound in the Wasserstein distance $d_{W}$, defined in (17), between $W_{n}$ and its limiting Beta distribution in (1). We show in Remark 3.1 that the rate of Theorem 1.1 cannot be improved. Let $x\wedge y$ and $x\vee y$ denote the minimum and maximum of two real numbers $x$ and $y$, respectively.

For $\alpha\geq 1$ and $\beta\geq 1$ let $S_{n}$ be the number of white balls in $n$ draws from a Pólya-Eggenberger urn that initially contains $\alpha$ white and $\beta$ black balls. Then with $W_{n}=S_{n}/n$ and $Z\sim{\cal B}(\alpha/m,\beta/m)$,

where $b_{0}=b_{0}(\alpha/m,\beta/m)$ and $b_{1}=b_{1}(\alpha/m,\beta/m)$ are given in Lemma 3.4.

Connections between Theorem 1.1 and the work of Döbler are spelled out in Remark 3.2.

The ${\cal B}(1/2,1/2)$ distribution, also known as the Arcsine law, describes the asymptotic distribution of many quantities that arise naturally in the study of the simple symmetric random walk $T_{n}=X_{1}+\cdots X_{n}$, where $X_{1},\ldots,X_{n}$ are independent variables taking the values 1 and $-1$ with probability $1/2$. For instance, let $L_{2n}$ be the random variable

giving the last return time to zero up to time $2n$. Then, see ,

where $P(T_{2m}=0)=2^{-2m}{2m\choose m}$, the probability that the walk returns to zero at time $2m$.

In the limit, $(2n)^{-1}L_{2n}\rightarrow Z$ in probability, where $Z$ has the Arcsine distribution. It is often noted that this limiting result is somewhat counter intuitive in that the Arcsine density has greatest mass near the endpoints, and least mass in the center of the unit interval, whereas in a fair coin tossing game one might assume that players are more likely to spend equal time in the lead. Perhaps at least as remarkable is the fact that the number $U_{2n}$ of segments of the walk that lie above the $x$ axis, and $R_{2n}$, the first time the walk visits the terminal point $S_{2n}$, are all equal in distribution to ${\cal L}(L_{2n})$; see for a nice exposition.

Using the methods presented here that were available in a preprint of this article, Döbler presents a Wasserstein bound of order $O(1/n)$ without explicit constants between the distribution of $(2n)^{-1}U_{2n}$ and the limiting Arcsine. In Section 4, essentially by applying the bounds in Lemma 3.4, we are able to attach concrete constants to the result of , as well as show the rate of the bound is optimal.

Let $L_{2n}$ be the last return time to zero of a simple symmetric random walk of length of length $2n$ and let $Z$ have the Arcsine distribution. Then

The same bound holds with $L_{2n}$ replaced by $U_{2n}$ or $R_{2n}$. The $O(1/n)$ rate of the bound cannot be improved.

Beginning with the introduction by Stein of a ‘characterizing equation’ method for developing bounds in normal approximation, to date the method has been successfully applied to a large number of the classical distributions, including the Poisson , , Multinomial , Gamma ,, Geometric , Negative Binomial and Exponential , , , as well as to non classical distributions such as the PRR family of also based on Pólya type urn models. Here we further extend the range of Stein’s method by including the Beta distribution, focusing on its role as the limiting law of the fraction of white balls drawn from the Pólya-Eggenberger urn.

The application of Stein’s method here differs from the way it is usually applied in that we focus on the approximation of particular distributions whose exact forms are known, rather than develop a bound that applies to a class of complex distributions obtained by, say, summing random variables that obey weak moment and dependence conditions, as in the case of the central limit theorem. And indeed, though explicit formulas exist for the distributions we study, the need for their approximation arises regardless, as is the case also for, say, the ubiquitous use of the normal approximation for the binomial.

Urn models of the classical type, and generalizations including drawing multiple balls or starting new urns, have received considerable attention recently; see for example , and . Interest has partly been sparked by the ability of urn models to exhibit power-law limiting behaviour, which in turn has been a focus of network analysis, see for example and . Connections between urn models and binary search trees are clearly explained in . In particular, let $m=1$ and consider the initial state of the Pólya-Eggenberger urn as a rooted binary tree having $\alpha$ white and $\beta$ black leaves, or external nodes. At every time step one external node is chosen, uniformly, to duplicate, yielding a pair of leaves of the same colour. That is, the chosen external node becomes an internal node while two external nodes of the chosen colour are added. The rule for adding an additional white leaf to the tree at time $n$ clearly is the same as the rule for adding an additional white ball to the Pólya-Eggenberger urn for the case $m=1$, and hence the number of white leaves of the tree and white balls in the urn have same distribution. Many variations and extensions on this theme are possible. Another line of interest comes from edge reinforced random walks, because an infinite system of independent Pólya-Eggenberger urns can be used to represent edge reinforced random walks on trees, see .

Characterizing equations and generators

Stein’s method for distributional approximation is based on a characterization of the target approximating distribution. For the seminal normal case considered in , it was shown that a variable $Z$ has the standard normal distribution if and only if

To obtain our result, we compute the distance between the distribution of the fraction of white balls drawn from the Pólya-Eggenberger Urn and the Beta by comparing the operators that characterize them. Our approach in characterizing the urn distribution stems from what is known as the density method; see for instance, , or Section 13.1 of . In particular, recognizing the $-w$ in (5) as the ratio of $\phi^{\prime}(w)/\phi(w)$ where $\phi(w)$ is the standard normal density, one hopes to replace the term $-w$ by the ratio $p^{\prime}(w)/p(w)$ when developing the Stein equation to handle the distribution with density $p(w)$, and to apply similar reasoning when the distribution under study is discrete. Use of the density method in the discrete case, followed by the application of a judiciously chosen transformation, leads to the characterization of the Pólya-Eggenberger Urn distribution given in Lemma 2.1.

Another approach to construct characterizing equations is known as the generator method. A number of years following the publication of , the relationship between the characterizing equation (4) and the generator of the Ornstein-Uhlenbeck process

of which the normal is the unique stationary measure, was recognized in , where it was noted that that in some generality the process semi-group may be used to solve the Stein equation (5). Given this connection between Stein characterizations and generators it is natural to consider a stochastic process which has the given target as its stationary distribution when extending Stein’s method to handle a new distribution.

Regarding the use of this ‘generator’ method for extending the scope of Stein’s method to the Beta distribution, we recall that the Fisher Wright model from genetics, originating in the work in , and , is a stochastic process used to model genetic drift in a population and has generator given by

Lemma 2.1 provides a characterizing equation for the Pólya urn distribution that is parallel to equation (6). Taking differences then allows us to estimate the expectation of the right hand side of (6) when $w$ is replaced by $W_{n}$ by exploiting the similarity of the two characterizing operators; a similar argument can be found in and for stationary distributions of birth-death chains. The results most closely related to the present work is , and its connections to the present manuscript are discussed in Remark 3.2

Let $p$ be the probability mass function of the number $S_{n}^{\alpha,\beta,m}$ of white balls drawn from the Pólya-Eggenberger urn by time $n$. Then a random variable $S$ has probability mass function $p$ if and only if for all functions $f\in{\cal{F}}(p)$

We prove Lemma 2.1 by applying a general technique for constructing equations such as (8) from discrete probability mass functions which is of independent interest, see . We begin with Proposition 2.1 below, a discrete version of the density approach to the Stein equation.

Let $Z$ have probability mass function $p$ with support the integer interval $I$, and let $\psi(k)$ be given by (7) for $k\in I$. If a random variable $X$ with support $I$ has mass function $p$ then for all $f\in{\cal{F}}(p)$,

The statement in Proposition 2.1 is equivalent to Theorem 1.1 given in under a different assumption, namely that equality (9) holds with $g$ replacing $f$ for all functions for which $\sum_{k\in I}\Delta(g(k)p(k))=0$. We note that their set-up would translate to test functions $f(k)=g(k+1)$.

If $p(b+1)=0$ and $f(a-1)=0$ then we obtain

Since $p(b+1)=0$ and $p(a)f(a-1)\rightarrow 0$ as $a\rightarrow-\infty$, we obtain that

Hence, if $I$ is $[a,b]$ or $[a,\infty)$ we obtain that for all $j\in I$,

Summing over $j\in I$ yields $P(X=a)=P(Z=a)$, whence $P(X=j)=P(Z=j)$ for all $j\in I$. Similarly one may handle the remaining case where $I=(-\infty,b]$. $\Box$

Given a characterization produced by Proposition 2.1, the following corollary produces varieties of characterizations for the same distribution, each one corresponding to a choice of a function $c$ possessing certain mild properties.

a nonstandard version of a characterization of the Poisson. An extension of Corollary 2.1 to the case of infinite support produces the usual characterization by the choice $c(k)=k+1$ and the substitution $g(k)=f(k-1)$. Naturally, additional characterizations are produced when using different choices of $c$.

where $(x)_{0}=1$ and otherwise $(x)_{k}=x(x+1)\cdots(x+k-1)$ is the rising factorial. The distribution (14) is also known as the beta-binomial and the negative hypergeometric distribution, see . We now have the ingredients to prove Lemma 2.1.

Proof of Lemma 2.1: Taking differences in (14) for $k=0,\ldots,n-1$ yields

Hence with $\psi(k)=\Delta p_{k}/p_{k}$ as in (7) we obtain for $k=0,\ldots,n-1$

In applying Corollary 2.1, as $\psi(n)=-1$ we may take the value $c(n)$ arbitrarily, see Remark 2.2. In particular, taking $c(k)=(k+1)(\beta/m+n-k-1)$ for all $k=0,\ldots,n-1$ and $c(n)=n$ we obtain (8). $\Box$

The next lemma is instrumental in calculating the higher moments of $S_{n}^{\alpha,\beta,m}$. We let $[x]_{0}=1$, and otherwise set $[x]_{k}=x(x-1)\cdots(x-k+1)$, the falling factorial.

For all nonnegative integers $n,a$ and $b$, we have

Proof: First we note that both sides of (15) are zero when $a+b\geq n+1$. This is clear for the right hand side, as the falling factorial $[n]_{a+b}$ is zero. For the left hand side, if $S_{n}\leq a-1$ then $[S_{n}]_{a}=0$. On the other hand, if $S_{n}\geq a$ then $b-1\geq n-a\geq n-S_{n}$, in which case $[n-S_{n}]_{b}$ is zero.

Now assume $n\geq a+b$. For any $k=0,1,\ldots,n$ we have

Summing over $k=0,1,\ldots,n$ and using that the support of $S_{r}$ is $\{0,\ldots,r\}$ yields (15). $\Box$

If $Z$ has the limiting beta distribution ${\cal B}(\alpha/m,\beta/m)$ with density (2), using (15) we obtain

that is, the scaled falling factorial moments of $S_{n}$ and the power moments of $Z$ differ only by factors of order $1/n$. This observation can be used to provide a proof of convergence in distribution of $W_{n}=S_{n}/n$ to $Z$ by the method of moments, but without a bound on the distributional distance.

Bounds for the Pólya-Eggenberger urn model

Theorem 1.1 provides an explicit bound in Wasserstein distance of order $O(1/n)$ between the distribution of $W_{n}$, the fraction of white balls drawn from the urn by time $n$, and the limiting Beta distribution. For approximating a discrete distribution by a continuous one the Wasserstein distance $d_{W}$ is a typical distance to use, see for example . For random variables $X$ and $Y$, this distance is given by

The function $h(x)=x(1-x){\bf 1}_{\{x\in\}}$ is in ${\rm Lip}(1)$, and applying (16) with $a=b=1$ we obtain that for all $\alpha\geq 1,\beta\geq 1$ and $m\geq 1$,

Thus the $1/n$ order of the bound in Theorem 1.1 cannot be improved.

Theorem 4.3 of provides a bound of order $1/n$ for the Beta approximation to the Pólya-Eggenberger urn for test functions with bounded first and second derivatives using an exchangeable pair coupling. The results in differ from ours in two significant ways. Firstly, the bound in Theorem 4.3 of is expressed in terms of two non-explicit constants $C_{1},C_{2}$ that are defined in Proposition 3.8 of . Lemma 3.4 below provides values of $C_{1}$. The lack of an explicit expression for $C_{2}$ in can be explained by the fact that the solution there is given in terms of ratios of functions which are related to incomplete Beta functions, for which a uniform bound would be difficult.

A more important difference between the present work and is that expressing the bound of the latter, presently given in terms of twice differentiable functions, in terms of a bound in a metric, say $d_{2}$, obtained from twice differentiable functions in the same way that Lipschitz functions yield the Wasserstein metric $d_{W}$, we have that $d_{2}\leq d_{W}$ with equality everywhere not holding. Hence Theorem 1.1 implies bounds in the $d_{2}$ metric, while the reverse does not hold.

In the following we set our test functions $h$ to be zero outside the unit interval $$. For$y>0$ set

and for a real valued function $g$ on $$we let$||g||=\sup_{w\in}|g(w)|$, the supremum norm of$g$. In the following we recall, with the help of Rademacher’s Theorem, that a function$h$is in${\rm Lip}(1)$ if and only it is absolutely continuous with respect to Lebesgue measure with an almost everywhere derivative bounded in absolute value by 1.

Lemma 3.1 below shows that for all $\{\alpha,\beta\}\subset(0,\infty)$ and functions $h$ for which the expectation ${\cal B}_{\alpha,\beta}h$ exists,

Proof of Theorem 1.1. For $h$ a given function in ${\rm Lip}(1)$, let $f=f_{\alpha/m,\beta/m}$ be the solution of the Stein equation (6) given in (18). Replacing $f(z)$ by $f(z/n)$ and dividing by $n$ in (8) results in

Applying this identity in the Stein equation (6), with $\alpha$ and $\beta$ replaced by $\alpha/m$ and $\beta/m$ respectively, and invoking Lemma 3.4 below to yield the existence and boundedness of $f^{\prime}$, we obtain

where, using Lemma 2.2 to calculate moments, we obtain

Writing the difference in (19) as an integral, we have

To handle $R_{2}$, using that the solution $f$ of the Stein equation equals for $x\not\in$ to obtain the first inequality,

For the first term in (21), substituting using the Stein equation (6) with $\alpha$ and $\beta$ replaced by $\alpha/m$ and $\beta/m$, respectively, we obtain

We bound the inner integrals separately. Firstly,

Next, recalling that $0\leq W_{n}\leq 1$ and noting that $|(\beta y-\alpha(1-y))|\leq\alpha\vee\beta$ for $0\leq y\leq 1$,

Collecting the bounds (20), (22), (23), (24) and (25) yields

The theorem now follows by invoking Lemma 3.4. $\Box$

For any $\{\alpha,\beta\}\subset(0,\infty)$ and real valued function $h$ on $$such that the expectation${\cal B}_{\alpha,\beta}h$of$h$exists, the function$f$ given by (18) is the unique bounded solution of (6).

Proof: It is straightforward to verify that $f$ as given in (18) is a solution of (6). Writing the associated homogeneous equation as

we find that all solutions to (6) are given by

The claim follows since $g(w)$ is unbounded at the endpoints of the unit interval for all $c\not=0$, and Lemma 3.4 below demonstrates that $f(w)$ is bounded. $\Box$

Since the expectation of $h(Z)-{\cal B}h$ is zero when $Z\sim{\cal B}(\alpha,\beta)$, we may also write

From Proposition 3.8 in we quote the following result.

The solution $f$, given in (18), of (6) for $h$ a Lipschitz function on $$ satisfies

The cases in the bounds of Lemma 3.4 reflect the behaviour of the function

as described in Lemma 3.3. In the following we will use the terms decreasing and increasing in the non-strict manner, for example, a constant function is both increasing and decreasing. Let

For $\{\alpha,\beta\}\subset(-1,\infty)$, the function $g:\rightarrow[0,\infty)$ given in (27) has the following behaviour.

Proof: Clearly when $\alpha=1$ and $\beta=1$ the function $g(w)$ is constant. Otherwise, taking derivative in (27) yields

The expression is non-negative if and only if

When $\alpha\geq 1$ and $\beta\leq 1$ inequality (28) is always satisfied. Similarly (28) holds with the non-strict inequality reversed when $\alpha\leq 1$ and $\beta\geq 1$. The remaining two cases $\alpha<1,\beta<1$ and $\alpha>1,\beta>1$ follow by solving the inequality. $\Box$

Our next result bounds the magnitude of the derivative of the solution $f$ in terms of $h$.

For $\{\alpha,\beta\}\subset(0,\infty)$ let $f=f_{\alpha,\beta,h}$ be the solution to (6) given by (18) for an absolutely continuous function $h$. Then

where $b_{0}=b_{0}(\alpha,\beta)$ and $b_{1}=b_{1}(\alpha,\beta)$ are given by

Proof: By replacing $h$ by $h-{\cal B}_{\alpha,\beta}h$ we may assume ${\cal B}_{\alpha,\beta}h=0$. Rewriting the Stein equation (6) yields

so to show (29) it suffices to demonstrate that for all $w\in$

Using (18) and integration by parts we obtain

From Lemma 3.2 we immediately have the bounds

When $u^{\alpha}(1-u)^{\beta-2}$ is increasing on $[0,x_{*}]$ then for $w\in[0,x_{*}]$ we can bound this expression by

and now using the first inequality in (32), we obtain

and if $u^{\alpha-2}(1-u)^{\beta}$ is decreasing on $[x_{*},1]$ then for $w\in[x_{*},1]$ we can bound this expression by

Now using the second inequality in (32), we obtain

In view of Lemma 3.3 we distinguish four cases.

Case 1. $\alpha\leq 2,\beta\leq 2$. By Lemma 3.3, $u^{\alpha}(1-u)^{\beta-2}$ is increasing and $u^{\alpha-2}(1-u)^{\beta}$ is decreasing. Setting $x_{*}=1/2$, by (33) and (34) we obtain

Case 2. $\alpha>2,\beta\leq 2$. In this case, from Lemma 3.3, $u^{\alpha}(1-u)^{\beta-2}$ is increasing, and $u^{\alpha-2}(1-u)^{\beta}$ is decreasing on $[x_{\alpha-1,\beta+1},1]$. Setting $x_{*}=x_{\alpha-1,\beta+1}$ and noting that $x_{\alpha,\beta}+x_{\beta,\alpha}=1$, by (33) and (34) we obtain

and bounding $(\alpha+\beta-2)/(\alpha+\beta)$ by 1 gives the assertion.

Case 3. $\alpha\leq 2,\beta>2$. In this case, from Lemma 3.3, $u^{\alpha}(1-u)^{\beta-2}$ is increasing on $[0,x_{\alpha+1,\beta-1}]$, and $u^{\alpha-2}(1-u)^{\beta}$ is decreasing. Setting $x_{*}=x_{\alpha+1,\beta-1}$, by (33) and (34) we obtain

Case 4. $\alpha>2,\beta>2$. In this case, from Lemma 3.3, $u^{\alpha}(1-u)^{\beta-2}$ is increasing on $[0,x_{\alpha+1,\beta-1}]$, and $u^{\alpha-2}(1-u)^{\beta}$ is decreasing on $[x_{\alpha-1,\beta+1},1]$. Noting that

setting $x_{*}=x_{\alpha,\beta}$, by (33) and (34) we obtain

For the final inequality in (29), with $p(y;\alpha,\beta)$ denoting the ${\cal B}(\alpha,\beta)$ density in (2), we have

We rely on for the following argument, noting that in no explicit bound is obtained.

Proof of Theorem 1.2. Let $p$ be the mass function of $L_{2n}$ given by (3), and let $Z$ have the Arcsine distribution. Applying Proposition 2.1 for $p$, followed by Corollary 2.1 with the choice

arrives at the version of Lemma 2.1, showing that $W_{n}=(2n)^{-1}L_{2n}$ is so distributed if and only if

Now following steps as those in Theorem 1.1 for the Pólya urn, collecting the estimates from the proof of Theorem 3.1 of shows that if $f$ is the solution (18) to (6) with $\alpha=\beta=1/2$, so that $b_{0}=2$ and $b_{1}=6$, then for all differentiable functions $h$ one has

Applying the bounds of Lemma 3.4 as well as Lemma 3.2 yields the bound in Theorem 1.2.

Applying (35) with $f(w)$ replaced by $g(w)=1$ and $g(w)=w$ yields

Hence the $O(1/n)$ rate cannot be improved. $\Box$

Acknowledgements. We would like to thank the Keble Advanced Studies Centre, Oxford, for support, and an anonymous referee for helpful comments which lead to an improvement of the paper.

References