TAOCP 3.4.1 Exercise 27

Let a subroutine $\operatorname{Bin}(m,\tfrac12)$ be available; it returns a random variable having the binomial distribution $(m,\tfrac12)$.

taocpmathematicsalgorithmsvolume-2medium

Section 3.4.1: Numerical Distributions

Exercise 27. [22] (J. H. Ahrens.) On most binary computers there is an efficient way to count the number of 1s in a binary word (see Section 7.1.3). Hence there is a nice way to obtain the binomial distribution $(t, p)$ when $p = \frac{1}{2}$, simply by generating $t$ random bits and counting the number of 1s.

Design an algorithm that produces the binomial distribution $(t, p)$ for arbitrary $p$, using only a subroutine for the special case $p = \frac{1}{2}$ as a source of random data. [Hint: Simulate a process that first looks at the most significant bits of $t$ uniform deviates, then at the second bit of those deviates whose leading bit is not sufficient to determine whether or not their value is $< p$, etc.]

Verified: yes
Solve time: 1m46s

Solution

Let a subroutine $\operatorname{Bin}(m,\tfrac12)$ be available; it returns a random variable having the binomial distribution $(m,\tfrac12)$. Equivalently, it supplies the number of 1s among $m$ independent random bits.

Write the binary expansion of $p$ as

$$ p=.b_1b_2b_3\cdots=\sum_{j\ge1}\frac{b_j}{2^j}, \qquad b_j\in{0,1}. $$

Consider a single uniform deviate

$$ U=.u_1u_2u_3\cdots . $$

The event $U<p$ is determined by comparing the binary expansions of $U$ and $p$ from left to right. Let

$$ A_j={u_1=b_1,\ldots,u_{j-1}=b_{j-1}}. $$

Conditional on $A_j$, the first $j-1$ bits do not determine whether $U<p$; the decision is made at bit $j$.

If $b_j=1$, then among the numbers satisfying $A_j$,

$$ u_j=0 \implies U<p, $$

while $u_j=1$ leaves the comparison undecided and requires examination of later bits.

If $b_j=0$, then among the numbers satisfying $A_j$,

$$ u_j=1 \implies U>p, $$

while $u_j=0$ leaves the comparison undecided.

Now suppose that $U_1,\ldots,U_t$ are independent uniform deviates. Let $R_j$ denote the number whose first $j-1$ bits agree with those of $p$; initially

$$ R_1=t. $$

At stage $j$, generate

$$ H_j\sim \operatorname{Bin}(R_j,\tfrac12). $$

Interpret $H_j$ as the number among the $R_j$ unresolved deviates whose $j$th bit equals $1$; hence $R_j-H_j$ have $j$th bit equal to $0$.

If $b_j=1$, then $R_j-H_j$ of the unresolved deviates are now known to satisfy $U_i<p$, while the $H_j$ deviates remain unresolved. Therefore set

$$ X\leftarrow X+(R_j-H_j),\qquad R_{j+1}\leftarrow H_j. $$

If $b_j=0$, then the $H_j$ deviates are known to satisfy $U_i>p$, while the remaining $R_j-H_j$ stay unresolved. Therefore set

$$ R_{j+1}\leftarrow R_j-H_j. $$

The process continues through the binary digits of $p$. Whenever $R_j=0$, all deviates have been classified and the algorithm terminates.

The variable $X$ counts the number of indices $i$ for which $U_i<p$.

To prove correctness, let

$$ I_i= \begin{cases} 1,&U_i<p,\ 0,&U_i\ge p. \end{cases} $$

Since $U_i$ is uniform on $[0,1)$,

$$ \Pr{I_i=1}=p. $$

The variables $I_1,\ldots,I_t$ are independent because the $U_i$ are independent. Hence

$$ X=\sum_{i=1}^t I_i $$

has the binomial distribution $(t,p)$.

The algorithm above merely determines the values of the indicators $I_i$ by revealing the bits of the uniforms sequentially. At stage $j$, the unresolved deviates are exactly those whose first $j-1$ bits coincide with the first $j-1$ bits of $p$; among them the next bits are independent Bernoulli variables with parameter $\tfrac12$. Therefore the number of next bits equal to $1$ is distributed as $\operatorname{Bin}(R_j,\tfrac12)$, exactly as generated by the available subroutine. Each stage preserves the correct conditional distribution of the unresolved deviates, and every deviate is eventually classified as either $U_i<p$ or $U_i>p$ (except on the null event $U_i=p$ when $p$ has a terminating binary expansion).

Consequently the final count $X$ is precisely the number of uniforms less than $p$, and therefore

$$ \Pr{X=k} = \binom{t}{k}p^k(1-p)^{t-k}, \qquad 0\le k\le t. $$

Thus the algorithm produces the binomial distribution $(t,p)$ using only calls to a subroutine for the special case $(m,\tfrac12)$. This completes the proof.

∎

Notes

A compact description is:

  1. Let $X\leftarrow0$, $R\leftarrow t$.
  2. For $j=1,2,\ldots$:
  • Generate $H\sim\operatorname{Bin}(R,\tfrac12)$.
  • If $b_j=1$, set $X\leftarrow X+(R-H)$ and $R\leftarrow H$.
  • If $b_j=0$, set $R\leftarrow R-H$.
  • Stop when $R=0$.

The expected value of $R_j$ is

$$ E(R_j)=t,2^{-(j-1)}, $$

so the number of stages required is typically of order $\log_2 t$.