TAOCP 3.4.1 Exercise 18

We are asked to generate a random integer $N$ such that \Pr\{N=n\} = n p^2 (1-p)^{\,n-1}, \qquad n \ge 0, with particular interest in the case where $p$ is small.

taocpmathematicsalgorithmsvolume-2math-medium

Section 3.4.1: Numerical Distributions

Exercise 18. [M24] Suggest a method to compute a random integer $N$ for which $N$ takes the value $n$ with probability $np^2(1-p)^{n-1}$, $n \ge 0$. (The case of particular interest is when $p$ is rather small.)

Verified: yes
Solve time: 4m33s

Exercise 3.4.1.18 [M24]

We are asked to generate a random integer $N$ such that

$$ \Pr{N=n} = n p^2 (1-p)^{,n-1}, \qquad n \ge 0, $$

with particular interest in the case where $p$ is small.

Step 1: Recognize the size-biased geometric distribution

Let $X$ be a geometric random variable with parameter $p$, taking values $n \ge 1$:

$$ \Pr{X=n} = p (1-p)^{,n-1}, \qquad n \ge 1. $$

Its expected value is

$$ \mathbb{E}[X] = \frac{1}{p}. $$

The target probability mass function can be written as

$$ \Pr{N=n} = n p^2 (1-p)^{,n-1} = p \cdot n p (1-p)^{,n-1} = \frac{n \Pr{X=n}}{\mathbb{E}[X]}, \qquad n \ge 1. $$

This shows that $N$ is the size-biased version of $X$: outcomes are weighted proportionally to $n$.

Thus, generating $N$ reduces to generating a size-biased geometric deviate.

Step 2: Construct the size-biased geometric distribution

Let $X$ be geometric as above. Define $N$ via the inverse transform method for size-biased distributions. The cumulative probabilities for $N$ are

$$ \Pr{N \le n} = \sum_{k=1}^{n} k p^2 (1-p)^{,k-1}. $$

To generate $N$, let $U \sim \text{Uniform}[0,1)$. We want $N$ satisfying

$$ \sum_{k=1}^{N-1} k p^2 (1-p)^{,k-1} \le U < \sum_{k=1}^{N} k p^2 (1-p)^{,k-1}. $$

This direct approach is inefficient for small $p$ because the sums can require many terms. A more efficient method uses two independent geometric deviates.

Step 3: Efficient algorithm using two independent geometric deviates

Let $X_1, X_2$ be independent geometric random variables with parameter $p$, taking values $n \ge 1$. Define

$$ N = X_1 + X_2 - 1. $$

We compute its distribution. Let $n \ge 1$. Then

$$ \begin{aligned} \Pr{N=n} &= \Pr{X_1 + X_2 - 1 = n} \ &= \Pr{X_1 + X_2 = n+1} \ &= \sum_{k=1}^{n} \Pr{X_1=k} \Pr{X_2=n+1-k} \ &= \sum_{k=1}^{n} p(1-p)^{,k-1} \cdot p(1-p)^{,n-k} \ &= p^2 \sum_{k=1}^{n} (1-p)^{,n-1} \ &= n p^2 (1-p)^{,n-1}. \end{aligned} $$

Thus $N$ defined as $X_1 + X_2 - 1$ has exactly the desired pmf. The case $n=0$ occurs with probability

$$ \Pr{N=0} = \Pr{X_1=1, X_2=1} = p^2 (1-p)^0 = p^2, $$

so the formula holds for $n \ge 0$.

Step 4: Algorithm implementation

To generate $N$:

  1. Generate two independent uniform random variables $U_1, U_2 \sim \text{Uniform}[0,1)$.
  2. Compute independent geometric deviates $X_i$ via inverse transform:

$$ X_i = \left\lfloor \frac{\ln(1-U_i)}{\ln(1-p)} \right\rfloor + 1, \quad i=1,2. $$

  1. Set

$$ N \leftarrow X_1 + X_2 - 1. $$

Then $N$ has the distribution

$$ \Pr{N=n} = n p^2 (1-p)^{,n-1}, \qquad n \ge 0. $$

This method is efficient for small $p$, requiring only two geometric deviates per sample.

Step 5: Summary

  • Recognize that the target law is a size-biased geometric distribution.
  • Use the fact that the sum of two independent geometric deviates minus one produces the size-biased pmf.
  • Generate two geometric deviates and combine them: $N = X_1 + X_2 - 1$.
  • This method is exact, simple, and efficient.

∎

This construction corrects the previous solution by properly realizing the size-biased geometric distribution, avoiding the invalid “minimum-of-two” method.