On the Average Runtime of an Open Source Binomial Random Variate Generation Algorithm

Abstract

The BTPE algorithm (Binomial, Triangle, Parallelogram, Exponential) of Kachitvichyanukul and Schmeiser is one of the faster and more widely utilized algorithms for generating binomial random variates. Cicirello’s open source Java library, $\rho\mu$, includes an implementation of BTPE as well as a variety of other random number related utilities. In this report, I explore the average case runtime of the BTPE algorithm when generating random values from binomial distribution $B(n,p)$. Beginning with Kachitvichyanukul and Schmeiser’s formula for the expected number of acceptance-rejection sampling iterations, I analyze the limit behavior as $n$ approaches infinity, and show that the average runtime of BTPE converges to a constant. I instrument the open source Java implementation from the $\rho\mu$ library to experimentally validate the analysis.

1 Introduction

In this report, I explore the average runtime behavior of binomial random variate generation within the open source Java library $\rho\mu$ [5]. The $\rho\mu$ library provides a variety of random number related enhancements beyond what is provided by the Java API itself. The core functionality of $\rho\mu$ is provided through a hierarchy of wrapper classes. This set of $\rho\mu$’s wrapper classes directly corresponds to the hierarchy of random number generator interfaces introduced in Java 17 [19]. In some cases, $\rho\mu$ overrides Java’s random number generators with faster algorithms, such as for random integers subject to a bound. In other cases, $\rho\mu$ adds functionality, such as additional distributions, i.e., the binomial, Cauchy, among others [5]. The initial motivation of the $\rho\mu$ library was to provide a source of enhanced and more efficient randomness for other libraries, such as JavaPermutationTools [3] and Chips-n-Salsa [4].

Among the randomness enhancements of the $\rho\mu$ library is support for generating binomial random variates, which are important to many applications [18, 17, 9, 22, 15, 2, 12, 21, 14]. There are many algorithms for generating binomial random variates [10, 11, 16, 1, 14, 13]. The $\rho\mu$ [5] library implements the BTPE Algorithm (Binomial, Triangle, Parallelogram, Exponential) [10], and falls back upon the inverse transform [10, 13] for cases that cannot be handled by BTPE.

The runtime of many algorithms for generating random values from binomial distribution $B(n,p)$ grows with some function of $n$ and $p$. For example, the runtime of the inverse transform method is $O(np)$ [10, 13]. One technique many binomial random variate algorithms utilize is acceptance–rejection sampling [8], and the iterations of such sampling is often what leads to such runtimes. BTPE does use acceptance–rejection sampling, but seems to require relatively few iterations. At the time of introduction, Kachitvichyanukul and Schmeiser determined the expected number of iterations as a function of $n$ and $p$ [10] (see Section 2.1). However, the aim of my present report is to provide simpler insight into the average case runtime than is provided by Kachitvichyanukul and Schmeiser’s formula.

In this report, I explore the limit of Kachitvichyanukul and Schmeiser’s formula as $n$ approaches infinity (Section 3), for two cases of $p$, including the minimum $p$ supported by BTPE (i.e., $p=\frac{10}{n}$), which maximizes the expected number of rejection sampling iterations, and $p=0.5$, which is when the expected number of rejection sampling iterations is minimized. Both cases lead to average runtimes that are constants in the limit for large $n$. In Section 4, I validate the findings experimentally using the $\rho\mu$ library, and utilizing my prior empirical study [6]. I discuss the findings and offer conclusions in Section 5, including the main result that the average runtime of BTPE is $\Theta(1)$. I now begin in Section 2 with some required background.

2 Preliminaries

2.1 Rejection Sampling Iterations of BTPE

The BTPE [10] algorithm separates the binomial distribution $B(n,p)$ into four parts, using triangular functions in the middle portions, and exponential functions in the tails, and it uses acceptance–rejection sampling [8]. For complete details of BTPE, which are beyond the scope of this paper, I refer the reader to the article that introduced it [10].

To generate a random value from binomial distribution $B(n,p)$, each acceptance–rejection iteration of BTPE generates two random values from $U(0,1)$, i.e., uniformly distributed over the interval $[0.0,1.0)$. When they introduced BTPE, Kachitvichyanukul and Schmeiser determined that the expected number of iterations, $E[I]$ of BTPE is [10]:

$$E[I]=p_{4}\binom{n}{M}r^{M}(1-r)^{n-M},$$ (1)

and since each iteration generates two random uniform values from $U(0,1)$, the expected number of uniform variates, $E[V]$, required by BTPE is thus:

$$E[V]=2p_{4}\binom{n}{M}r^{M}(1-r)^{n-M}.$$ (2)

Equations 1 and 2 involve $r$, $M$, and $p_{4}$, which in turn depend upon a few others, all of which are functions of $n$ and $p$. Kachitvichyanukul and Schmeiser [10] define these as follows:

$$r=\min(p,1-p),$$ (3)

$$q=1-r,$$ (4)

$$f_{M}=nr+r,$$ (5)

$$M=\lfloor f_{M}\rfloor,$$ (6)

$$p_{1}=\left\lfloor 2.195\sqrt{nrq}-4.6q\right\rfloor+0.5,$$ (7)

$$x_{M}=M+0.5,$$ (8)

$$x_{L}=x_{M}-p_{1},$$ (9)

$$x_{R}=x_{M}+p_{1},$$ (10)

$$c=0.134+\frac{20.5}{15.3+M},$$ (11)

$$a_{L}=\frac{f_{M}-x_{L}}{f_{M}-x_{L}r},$$ (12)

$$\lambda_{L}=a_{L}\left(1+\frac{a_{L}}{2}\right),$$ (13)

$$a_{R}=\frac{x_{R}-f_{M}}{x_{R}q},$$ (14)

$$\lambda_{R}=a_{R}\left(1+\frac{a_{R}}{2}\right),$$ (15)

$$p_{2}=p_{1}(1+2c),$$ (16)

$$p_{3}=p_{2}+\frac{c}{\lambda_{L}},$$ (17)

$$p_{4}=p_{3}+\frac{c}{\lambda_{R}}.$$ (18)

Also note that BTPE is only relevant when: $\min(p,1-p)\geq\frac{10}{n}$. Kachitvichyanukul and Schmeiser recommend using the inverse transform method [13] as a fallback for cases when BTPE is not relevant. This is what is done in the implementation in the $\rho\mu$ library [5].

The expected number of acceptance–rejection sampling iterations, expressed in Equation 2, is maximized when $\min(p,1-p)$ is minimal. Thus, the expected number of rejection sampling iterations is maximized for the cases $p=\frac{10}{n}$ and $p=\frac{n-10}{n}$. Due to symmetry, these two cases lead to the same expected number of rejection sampling iterations. The nearer $p$ is to $0.5$, the fewer iterations BTPE requires on average. Table 1 shows how the expected number of required uniform variates changes as $n$ increases for these two cases of $p$. It includes the minimum $n$ supported by BTPE (e.g., $n=20$), and then considers $n$ at increasing powers of two to show how the number of uniform variates used by BTPE to generate one binomial random variate varies as $n$ increases rapidly.

Table 1: Expected number of uniform variates to generate one binomial random variate as predicted by Equation 2 as $n$ increases.

	$E[V]$ from Equation 2
$n$	$p\in\left\{\frac{10}{n},\frac{n-10}{n}\right\}$	$p=0.5$
$20$	3.996	3.996
$2^{5}$	3.837	3.599
$2^{6}$	3.792	3.337
$2^{7}$	3.790	2.985
$2^{8}$	3.793	2.632
$2^{9}$	3.796	2.420
$2^{10}$	3.797	2.317
$2^{11}$	3.798	2.307
$2^{12}$	3.798	2.276
$2^{13}$	3.799	2.300
$2^{14}$	3.799	2.299
$2^{15}$	3.799	2.300
$2^{16}$	3.799	2.302
$2^{17}$	3.799	2.310
$2^{18}$	3.799	2.310
$2^{19}$	3.799	2.313
$2^{20}$	3.799	2.314

2.2 Stirling’s Formula

While computing the limiting behavior of BTPE, we will encounter some factorials. Stirling’s formula [7] will be useful in the analysis. Stirling’s formula is as follows:

$$n!=\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}e^{\alpha_{n}},$$ (19)

where

$$\frac{1}{12n+1}<\alpha_{n}<\frac{1}{12n}.$$ (20)

A consequence of using Stirling’s formula is that it will introduce terms involving $e^{\alpha_{n}}$, either generally for $n$ or in other cases for specific values of $n$. The following observations will be useful:

$$\lim_{n\to\infty}\alpha_{n}=0\>,$$ (21)

$$\lim_{n\to\infty}\alpha_{n/2}=0\>,$$ (22)

$$\lim_{n\to\infty}\alpha_{n-10}=0\>,$$ (23)

$$\lim_{n\to\infty}e^{\alpha_{n}}=1\>,$$ (24)

$$\lim_{n\to\infty}e^{\alpha_{n/2}}=1\>,$$ (25)

$$\lim_{n\to\infty}e^{\alpha_{n-10}}=1\>.$$ (26)

3 Limit Analysis

Let’s now proceed to compute the limit:

$$\lim_{n\to\infty}E[V]=\lim_{n\to\infty}2p_{4}\binom{n}{M}r^{M}(1-r)^{n-M},$$ (27)

for two cases, when $p=\frac{10}{n}$ and when $p=0.5$ in Sections 3.1 and 3.2, respectively.

3.1 Case 1: Maximum Rejection Sampling Iterations

Consider the limit from Equation 27 when $p=\frac{10}{n}$, which is the lowest $p$ (in terms of $n$) supported by BTPE. Begin by computing the values of the various constants from Section 2.1 in terms of $n$ and $p=\frac{10}{n}$ as follows:

$$r=\min(p,1-p)=\frac{10}{n}\>,$$ (28)

$$q=1-r=\frac{n-10}{n}\>,$$ (29)

$$f_{M}=nr+r=n\left(\frac{10}{n}\right)+\frac{10}{n}=10+\frac{10}{n}\>,$$ (30)

$$M=\lfloor f_{M}\rfloor=\left\lfloor 10+\frac{10}{n}\right\rfloor=10\>,$$ (31)

$$c=0.134+\frac{20.5}{15.3+M}=0.134+\frac{20.5}{15.3+10}\approx 0.944\>,$$ (32)

$$x_{M}=M+0.5=10+0.5=10.5\>.$$ (33)

From Equation (20), find the following, which we need later:

$$\frac{1}{121}<\alpha_{10}<\frac{1}{120}\>,$$ (34)

$$\alpha_{10}\approx 0.0083\>,$$ (35)

$$e^{\alpha_{10}}\approx 1.008\>.$$ (36)

We will also need the following limit:

$$\begin{split}\lim_{n\to\infty}p_{1}&=\lim_{n\to\infty}\left\lfloor 2.195\sqrt{nrq}-4.6q\right\rfloor+0.5\\ &=\lim_{n\to\infty}\left\lfloor 2.195\sqrt{n\left(\frac{10}{n}\right)\left(\frac{n-10}{n}\right)}-4.6\left(\frac{n-10}{n}\right)\right\rfloor+0.5\\ &=\lim_{n\to\infty}\left\lfloor 2.195\sqrt{10\left(\frac{n-10}{n}\right)}-4.6\left(\frac{n-10}{n}\right)\right\rfloor+0.5\\ &=\left\lfloor 2.195\sqrt{10}-4.6\right\rfloor+0.5\\ &=2.5\>.\end{split}$$ (37)

We will also need the limits of $\lambda_{L}$ and $\lambda_{R}$ as $n$ approaches infinity, as follows:

$$\begin{split}\lim_{n\to\infty}\lambda_{L}&=\lim_{n\to\infty}a_{L}\left(1+\frac{a_{L}}{2}\right)\\ &=\lim_{n\to\infty}\left(\frac{f_{M}-x_{L}}{f_{M}-x_{L}r}\right)\left(1+\frac{f_{M}-x_{L}}{2(f_{M}-x_{L}r)}\right)\\ &=\lim_{n\to\infty}\left(\frac{f_{M}-x_{M}+p_{1}}{f_{M}-r(x_{M}-p_{1})}\right)\left(1+\frac{f_{M}-x_{M}+p_{1}}{2(f_{M}-r(x_{M}-p_{1}))}\right)\\ &=\lim_{n\to\infty}\left(\frac{f_{M}-8}{f_{M}-8r}\right)\left(1+\frac{f_{M}-8}{2(f_{M}-8r)}\right)\\ &=\lim_{n\to\infty}\left(\frac{\frac{10n+10}{n}-8}{\frac{10n+10}{n}-\frac{80}{n}}\right)\left(1+\frac{\frac{10n+10}{n}-8}{2(\frac{10n+10}{n}-\frac{80}{n})}\right)\\ &=\left(\frac{10-8}{10}\right)\left(1+\frac{10-8}{20}\right)\\ &=0.22\end{split}$$ (38)

and

$$\begin{split}\lim_{n\to\infty}\lambda_{R}&=\lim_{n\to\infty}a_{R}\left(1+\frac{a_{R}}{2}\right)\\ &=\lim_{n\to\infty}\left(\frac{x_{R}-f_{M}}{x_{R}q}\right)\left(1+\frac{x_{R}-f_{M}}{2x_{R}q}\right)\\ &=\lim_{n\to\infty}\left(\frac{x_{M}+p_{1}-f_{M}}{q(x_{M}+p_{1})}\right)\left(1+\frac{x_{M}+p_{1}-f_{M}}{2q(x_{M}+p_{1})}\right)\\ &=\lim_{n\to\infty}\left(\frac{13-f_{M}}{13q}\right)\left(1+\frac{13-f_{M}}{26q}\right)\\ &=\lim_{n\to\infty}\left(\frac{13-\frac{10n+10}{n}}{\frac{13n-130}{n}}\right)\left(1+\frac{13-\frac{10n+10}{n}}{\frac{26n-260}{n}}\right)\\ &=\left(\frac{13-10}{13}\right)\left(1+\frac{13-10}{26}\right)\\ &\approx 0.257.\end{split}$$ (39)

We finally can compute the limit of Equation (27):

$$\begin{split}\lim_{n\to\infty}E[V]&=\lim_{n\to\infty}2p_{4}\binom{n}{M}r^{M}(1-r)^{n-M}\\ &=\lim_{n\to\infty}2p_{4}\binom{n}{10}\left(\frac{10}{n}\right)^{10}\left(\frac{n-10}{n}\right)^{n-10}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{n!}{10!(n-10)!}\right)\left(\frac{10}{n}\right)^{10}\left(\frac{n-10}{n}\right)^{n-10}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{\sqrt{\frac{n}{(n-10)}}\,n^{n}\,e^{\alpha_{n}}}{\sqrt{20\pi}\,10^{10}(n-10)^{n-10}\,e^{\alpha_{n-10}}\,e^{\alpha_{10}}}\right)\left(\frac{10}{n}\right)^{10}\left(\frac{n-10}{n}\right)^{n-10}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{\sqrt{\frac{n}{(n-10)}}}{\sqrt{20\pi}\,e^{\alpha_{10}}}\right)\approx\lim_{n\to\infty}2p_{4}\left(\frac{\sqrt{\frac{n}{(n-10)}}}{1.008\sqrt{20\pi}}\right)\approx\lim_{n\to\infty}0.2503p_{4}\sqrt{\frac{n}{(n-10)}}\\ &\approx\lim_{n\to\infty}0.2503p_{4}\\ &\approx\lim_{n\to\infty}0.2503\left(p_{1}(1+2c)+\frac{c}{\lambda_{L}}+\frac{c}{\lambda_{R}}\right)\\ &\approx\lim_{n\to\infty}0.2503\left(2.888\,p_{1}+\frac{0.944}{\lambda_{L}}+\frac{0.944}{\lambda_{R}}\right)\\ &\approx\lim_{n\to\infty}0.2503\left(7.22+\frac{0.944}{\lambda_{L}}+\frac{0.944}{\lambda_{R}}\right)\\ &\approx 0.2503\left(7.22+\frac{0.944}{0.22}+\frac{0.944}{0.257}\right)\approx 0.2503\left(7.22+4.291+3.673\right)\\ &\approx 3.801.\end{split}$$ (40)

Thus, the expected number of uniform random variates required by BTPE to generate one binomial random variate converges to approximately 3.801 as $n$ grows large (approximately 1.90 rejection sampling iterations on average) for the case of minimum supported $p$. Observe in Table 1 that BTPE approaches this limit case rapidly.

3.2 Case 2: Central Case

Now consider the limit from Equation 27 for the case of $p=0.5$.

Let us begin by computing some of the constants from Section 2.1 for this case. When computing $M$, assume that $n$ is even, without loss of generality:

$$r=\min(p,1-p)=\frac{1}{2},$$ (41)

$$q=1-r=\frac{1}{2},$$ (42)

$$f_{M}=nr+r=\frac{n+1}{2},$$ (43)

$$M=\lfloor f_{M}\rfloor=\frac{n}{2}\>.$$ (44)

Now substitute these, and other Section 2.1 definitions, into Equation 27 and simplify:

$$\begin{split}\lim_{n\to\infty}E[V]&=\lim_{n\to\infty}2p_{4}\binom{n}{M}r^{M}(1-r)^{n-M}\\ &=\lim_{n\to\infty}2p_{4}\binom{n}{\frac{n}{2}}\left(\frac{1}{2}\right)^{\frac{n}{2}}\left(1-\frac{1}{2}\right)^{n-\frac{n}{2}}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{n!}{\left(\frac{n}{2}\right)!\left(\frac{n}{2}\right)!}\right)\left(\frac{1}{2}\right)^{n}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}e^{\alpha_{n}}}{\sqrt{\pi n}\left(\frac{n}{2e}\right)^{n/2}e^{\alpha_{n/2}}\sqrt{\pi n}\left(\frac{n}{2e}\right)^{n/2}e^{\alpha_{n/2}}}\right)\left(\frac{1}{2}\right)^{n}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}e^{\alpha_{n}}}{\pi n\left(\frac{n}{e}\right)^{n}\left(\frac{1}{2}\right)^{n}e^{2\alpha_{n/2}}}\right)\left(\frac{1}{2}\right)^{n}\\ &=\lim_{n\to\infty}2p_{4}\left(\frac{\sqrt{2/\pi}}{\sqrt{n}}\right)\left(\frac{e^{\alpha_{n}}}{e^{2\alpha_{n/2}}}\right)\\ &\approx\lim_{n\to\infty}1.596\left(\frac{p_{4}}{\sqrt{n}}\right)\\ &\approx\lim_{n\to\infty}1.596\left(\frac{1}{\sqrt{n}}\right)\left(p_{1}(1+2c)+\frac{c}{\lambda_{L}}+\frac{c}{\lambda_{R}}\right)\\ &\approx\lim_{n\to\infty}1.596\left(\frac{1}{\sqrt{n}}\right)\left(p_{1}(1+2c)+\frac{2c}{\lambda_{R}}\right)\\ &\approx\lim_{n\to\infty}1.596\left(\frac{p_{1}(1+2c)}{\sqrt{n}}+\frac{2c}{\sqrt{n}\lambda_{R}}\right).\end{split}$$ (45)

The next-to-last step above utilizes the fact that $\lambda_{R}=\lambda_{L}$ when $p=0.5$.

Before continuing to simplify Equation 45, compute the following limits. First:

$$\lim_{n\to\infty}c=\lim_{n\to\infty}0.134+\frac{20.5}{15.3+M}=\lim_{n\to\infty}0.134+\frac{20.5}{15.3+\frac{n}{2}}=0.134.$$ (46)

Now, compute the limit $\lim_{n\to\infty}\frac{p_{1}}{\sqrt{n}}$. Since $p_{1}$ involves a floor, compute bounds on this limit. Although we will see that the upper and lower bounds are equal, and thus via the squeeze theorem is also equal to our target limit.

$$\begin{split}\lim_{n\to\infty}\frac{\left\lfloor 2.195\sqrt{nrq}-4.6q\right\rfloor+0.5}{\sqrt{n}}\>&\leq\lim_{n\to\infty}\frac{p_{1}}{\sqrt{n}}\leq\>\lim_{n\to\infty}\frac{\left\lfloor 2.195\sqrt{nrq}-4.6q\right\rfloor+0.5}{\sqrt{n}}\\ \lim_{n\to\infty}\frac{\left\lfloor 2.195\sqrt{n/4}-2.3\right\rfloor+0.5}{\sqrt{n}}\>&\leq\lim_{n\to\infty}\frac{p_{1}}{\sqrt{n}}\leq\>\lim_{n\to\infty}\frac{\left\lfloor 2.195\sqrt{n/4}-2.3\right\rfloor+0.5}{\sqrt{n}}\\ \lim_{n\to\infty}\frac{\left\lfloor 1.0975\sqrt{n}-2.3\right\rfloor+0.5}{\sqrt{n}}\>&\leq\lim_{n\to\infty}\frac{p_{1}}{\sqrt{n}}\leq\>\lim_{n\to\infty}\frac{\left\lfloor 1.0975\sqrt{n}-2.3\right\rfloor+0.5}{\sqrt{n}}\\ \lim_{n\to\infty}\frac{1.0975\sqrt{n}-2.8}{\sqrt{n}}\>&\leq\lim_{n\to\infty}\frac{p_{1}}{\sqrt{n}}\leq\>\lim_{n\to\infty}\frac{1.0975\sqrt{n}-1.8}{\sqrt{n}}\\ 1.0975\>&\leq\lim_{n\to\infty}\frac{p_{1}}{\sqrt{n}}\leq\>1.0975\end{split}$$ (47)

The next limit that we need is $\lim_{n\to\infty}\frac{2}{\sqrt{n}\lambda_{R}}$ as follows:

$$\begin{split}\lim_{n\to\infty}\frac{2}{\sqrt{n}\lambda_{R}}&=\lim_{n\to\infty}\frac{2}{\sqrt{n}\left(a_{R}\left(1+\frac{a_{R}}{2}\right)\right)}\\ &=\lim_{n\to\infty}\frac{2x_{R}q}{\sqrt{n}(x_{R}-f_{M})\left(1+\frac{x_{R}-f_{M}}{2x_{R}q}\right)}\\ &=\lim_{n\to\infty}\frac{x_{R}}{\sqrt{n}(x_{R}-f_{M})\left(1+\frac{x_{R}-f_{M}}{x_{R}}\right)}\end{split}$$ (48)

Now consider the following:

$$x_{R}-f_{M}=x_{M}+p_{1}-\frac{n+1}{2}=M+\frac{1}{2}+p_{1}-\frac{n+1}{2}=\frac{n}{2}+\frac{1}{2}+p_{1}-\frac{n+1}{2}=p_{1}.$$ (49)

Now continue simplifying Equation 48, keeping Equation 47 in mind as well as that $p_{1}\in\Theta(\sqrt{n})$ as follows:

$$\begin{split}\lim_{n\to\infty}\frac{2}{\sqrt{n}\lambda_{R}}&=\lim_{n\to\infty}\frac{x_{R}}{\sqrt{n}(x_{R}-f_{M})\left(1+\frac{x_{R}-f_{M}}{x_{R}}\right)}\\ &=\lim_{n\to\infty}\frac{x_{M}+p_{1}}{p_{1}\sqrt{n}\left(1+\frac{p_{1}}{x_{M}+p_{1}}\right)}\\ &=\lim_{n\to\infty}\frac{\frac{n+1}{2}+p_{1}}{p_{1}\sqrt{n}\left(1+\frac{p_{1}}{\frac{n+1}{2}+p_{1}}\right)}\\ &=\lim_{n\to\infty}\frac{n+1+2p_{1}}{2p_{1}\sqrt{n}\left(1+\frac{2p_{1}}{n+1+2p_{1}}\right)}\\ &=\lim_{n\to\infty}\frac{n+1+2p_{1}}{2p_{1}\sqrt{n}}\\ &=\lim_{n\to\infty}\frac{\sqrt{n}+\frac{1}{\sqrt{n}}+2\frac{p_{1}}{\sqrt{n}}}{2p_{1}}\\ &=\lim_{n\to\infty}\frac{2.195+\sqrt{n}}{2p_{1}}\\ &=\lim_{n\to\infty}\frac{\sqrt{n}}{2p_{1}}=\frac{1}{2(1.0975)}=\frac{1}{2.195}\approx 0.456\>.\end{split}$$ (50)

Now continue simplifying Equation 45, utilizing Equations 46, 47, and 50 as follows:

$$\begin{split}\lim_{n\to\infty}E[V]&\approx\lim_{n\to\infty}1.596\left(\frac{p_{1}(1+2c)}{\sqrt{n}}+\frac{2c}{\sqrt{n}\lambda_{R}}\right)\\ &\approx\lim_{n\to\infty}1.596\left(1.0975(1+2c)+0.456c\right)\\ &\approx 1.596\left(1.0975\cdot 1.268+0.456\cdot 0.134\right)\\ &\approx 2.319.\end{split}$$ (51)

Thus, when $p=0.5$ the expected number of uniform random variates required by BTPE to generate one binomial random variate converges to approximately 2.319 as $n$ grows large (approximately 1.16 rejection sampling iterations on average). Observe in Table 1, that this limit case is reached quickly.

4 Experimental Validation

In earlier work, I used the open source Java library $\rho\mu$ to empirically explore the behavior of my BTPE implementation [6]. I wrapped an instance of Java’s SplittableRandom class, which implements the splitmix [20] pseudorandom number generator (PRNG). The purpose was to instrument the PRNG in order to count the number of uniform random variates, $U(0,1)$, generated (i.e., calls to the nextDouble() method) while generating a binomial random variate. I then supply this wrapped PRNG instance as the source of randomness for $\rho\mu$’s implementation of BTPE.

That prior paper [6] considered many combinations of $n$ and $p$. Here I focus on the same values of $n\in\{2^{5},2^{6},\ldots,2^{20}\}$. But, I only focus on two specific cases of $p$, the two cases utilized earlier in the paper: $p=\frac{10}{n}$ from Section 3.1 and $p=0.5$ from Section 3.2. For each binomial distribution $B(n,p)$, I generate 10,000 binomial random variates, and I compute the average number of uniform variates per binomial, with 95% confidence intervals. I compare the experimental mean to the prediction of the number of uniform variates from Equation 2, and test significance with a $t$-test. The experiments used OpenJDK 17 on a Windows 10 PC with a 3.4 GHz AMD A10-5700 CPU and 8 GB RAM. I used $\rho\mu$ 3.1.1. Both the source code for the experiments (https://github.com/cicirello/btpe-iterations), as well as for $\rho\mu$ (https://github.com/cicirello/rho-mu) is maintained on GitHub. The API documentation for the $\rho\mu$ library is available on the web: https://rho-mu.cicirello.org/.

Table 2: Average number of calls to $U(0,1)$ by $\rho\mu$’s BTPE implementation compared to the prediction of Equation 2. 95% confidence intervals are shown, as well as the $p$-values from $t$-tests.

	$B(n,\frac{10}{n})$			$B(n,0.5)$
$n$	predicted	mean	$p$-value	predicted	mean	$p$-value
$2^{5}$	3.837	$3.829\pm 0.052$	0.76	3.599	$3.582\pm 0.046$	0.46
$2^{6}$	3.792	$3.786\pm 0.050$	0.80	3.337	$3.334\pm 0.042$	0.89
$2^{7}$	3.790	$3.810\pm 0.051$	0.45	2.985	$3.004\pm 0.034$	0.27
$2^{8}$	3.793	$3.765\pm 0.051$	0.27	2.632	$2.637\pm 0.025$	0.66
$2^{9}$	3.796	$3.796\pm 0.051$	0.99	2.420	$2.412\pm 0.019$	0.41
$2^{10}$	3.797	$3.818\pm 0.051$	0.42	2.317	$2.304\pm 0.016$	0.11
$2^{11}$	3.798	$3.875\pm 0.052$	0.00	2.307	$2.308\pm 0.017$	0.85
$2^{12}$	3.798	$3.808\pm 0.051$	0.72	2.276	$2.276\pm 0.015$	0.98
$2^{13}$	3.799	$3.809\pm 0.051$	0.70	2.300	$2.288\pm 0.016$	0.13
$2^{14}$	3.799	$3.798\pm 0.051$	0.97	2.299	$2.304\pm 0.017$	0.55
$2^{15}$	3.799	$3.808\pm 0.052$	0.73	2.300	$2.286\pm 0.016$	0.08
$2^{16}$	3.799	$3.828\pm 0.052$	0.27	2.302	$2.306\pm 0.017$	0.65
$2^{17}$	3.799	$3.825\pm 0.052$	0.32	2.310	$2.301\pm 0.016$	0.29
$2^{18}$	3.799	$3.818\pm 0.051$	0.46	2.310	$2.313\pm 0.017$	0.73
$2^{19}$	3.799	$3.792\pm 0.050$	0.80	2.313	$2.314\pm 0.017$	0.91
$2^{20}$	3.799	$3.807\pm 0.051$	0.74	2.314	$2.309\pm 0.017$	0.58

Table 2 shows the results. The raw and processed data are available on GitHub (https://github.com/cicirello/btpe-iterations). The empirical results confirm the analytical prediction of Equation (2). Observe that there is no significant difference between the analytical prediction and the empirically computed means. T-test $p$-values are above 0.05 in almost all cases (well above in most cases). There is only a single case, $B(2^{11},10/2^{11})$, in Table 2 where a $t$-test $p$-value is less than 0.05, seemingly suggesting that the implementation behaves differently than theory predicts in this case. However, this case is explainable by random chance. There are 32 cases represented in Table 2, so this one case is approximately 3% of the cases tested. My more in depth empirical study explored a much larger number of cases experimentally, finding no significant difference between predicted and observed behavior in all but 4% of cases [6]. At significance level 0.05, we should expect false-positives in approximately 5% of cases.

5 Discussion and Conclusions

In this report, I analyzed the runtime behavior of the BTPE algorithm for binomial random variate generation, and observed the following concerning average performance:

• The average number of acceptance–rejection sampling iterations, $E[I]$, and average number of uniform random numbers generated while generating a binomial random variate, $E[V]$, are as follows:

  $$1\leq E[I]<2,$$ (52)

  $$2\leq E[V]<4.$$ (53)

  At least one iteration is necessary (the above lower bound), and the above upper bounds derive from the case of minimum supported $n=20$ for $p=0.5$, which is the case requiring the maximum number of rejection sampling iterations on average (Table 1).

• For large $n$ and the case of $p=\frac{10}{n}$, which maximizes the expected number of rejection sampling iterations, we find that the expected number of uniform random numbers required to generate one binomial random variate is a constant:

  $$\lim_{n\to\infty}E[V]\approx 3.801.$$ (54)

• For large $n$ and the case of $p=0.5$, the case that minimizes the expected number of rejection sampling iterations, we find that the expected number of uniform random numbers required to generate one binomial random variate is also a constant:

  $$\lim_{n\to\infty}E[V]\approx 2.319.$$ (55)

Thus, the average runtime of BTPE is $\Theta(1)$, while the average runtime of many of the alternative algorithms for generating binomial random variates is worse than a constant. For example, the average runtime of the inverse transform method is $O(np)$ [10, 13].

In this report, I also experimentally validated my limit analysis of Kachitvichyanukul and Schmeiser’s formula, using my implementation of BTPE in the $\rho\mu$ Java library [5]. This validates: (a) Kachitvichyanukul and Schmeiser’s formula for the expected number of rejection sampling iterations, (b) my analysis of that formula, and (c) that $\rho\mu$’s implementation of BTPE behaves as theory predicts. One limitation of this study is that it focuses on average case performance. It is possible that BTPE may exhibit many more rejection sampling iterations on occasion. The maximum number of uniform random variates generated while generating a single binomial during the course of my previous experimentation [6] was 38 (19 rejection sampling iterations). Such instances were very rare in that study, which generated approximately 2.88 million binomial random variates.