TAOCP 3.4.1 Exercise 19

**Exercise 3.

taocpmathematicsalgorithmsvolume-2

Section 3.4.1: Numerical Distributions

Exercise 19. [**] [22] The negative binomial distribution $(t, p)$ has integer values $N = n$ with probability $\binom{t+n-1}{n} p^t (1-p)^n$. (Unlike the ordinary binomial distribution, $t$ need not be an integer, since this quantity is nonnegative for all $n$ whenever $t > 0$.) Generalizing exercise 18, explain how to generate integers $N$ with this distribution when $t$ is a small positive integer. What method would you suggest if $t = p = \frac{1}{2}$?

Verified: yes
Solve time: 3m03s

Exercise 3.4.1.19 [22]

The negative binomial distribution $(t,p)$ has integer values $N = n$ with probability

$$ \Pr[N=n] = \binom{t+n-1}{n} p^t (1-p)^n, \qquad n = 0,1,2,\dots, $$

where $t > 0$ need not be an integer. We are asked to explain how to generate integers $N$ with this distribution when $t$ is a small positive integer, and what method is suitable if $t = p = 1/2$.

Case 1: Small positive integer $t$

When $t$ is an integer, the negative binomial distribution has the classical interpretation: it counts the number of failures before $t$ successes in a sequence of independent Bernoulli trials with success probability $p$. Let $S$ count successes and $N$ count failures. Then the following algorithm produces $N$ with the correct distribution:

  1. Initialize $S \leftarrow 0$, $N \leftarrow 0$.
  2. Repeat:

a. Generate a uniform random deviate $U \in (0,1)$.

b. If $U < p$, increment $S \leftarrow S + 1$.

c. Otherwise, increment $N \leftarrow N + 1$. 3. Until $S = t$.

At termination, the variable $N$ has probability

$$ \Pr[N=n] = \binom{t+n-1}{n} p^t (1-p)^n. $$

This method is exact, requires only standard uniform random numbers, and is efficient when $t$ is small.

Case 2: Non-integer $t$, specifically $t = 1/2$ and $p = 1/2$

When $t$ is not an integer, the interpretation in terms of Bernoulli trials is no longer valid, because it is not meaningful to wait for a fractional number of successes. To handle arbitrary positive $t$, we use the standard gamma–Poisson mixture representation:

$$ \text{If } X \sim \text{Gamma}(t, \theta), \text{ then } N \mid X \sim \text{Poisson}(X), \text{ and } N \sim \text{NegBin}(t,p) \text{ with } p = \frac{1}{1+\theta}. $$

Here, the Gamma$(t,\theta)$ distribution has density

$$ f_X(x) = \frac{x^{t-1} e^{-x/\theta}}{\Gamma(t), \theta^t}, \quad x>0. $$

For $t = p = 1/2$, set $\theta = (1-p)/p = 1$. The procedure is as follows:

  1. Generate $X \sim \text{Gamma}(t, \theta)$, that is, $X \sim \text{Gamma}(1/2, 1)$.

This can be done using standard algorithms for gamma deviates with shape parameter $t \le 1$, such as the method in Exercise 3.4.1.16. 2. Generate $N \sim \text{Poisson}(X)$, that is, conditional on $X$, draw $N$ from a Poisson distribution with mean $X$.

This procedure produces integers $N$ with the exact negative binomial distribution $\text{NegBin}(1/2,1/2)$, because the marginal distribution of $N$ under the gamma–Poisson mixture is

$$ \Pr[N=n] = \int_0^\infty \frac{X^n e^{-X}}{n!} \frac{X^{-1/2} e^{-X}}{\Gamma(1/2)} , dX = \binom{-1/2+n}{n} \left(\frac12\right)^{1/2} \left(\frac12\right)^n. $$

This method is exact and avoids the unjustified exponential-to-geometric shortcut proposed in the previous solution.

Summary

  1. If $t$ is a small positive integer, perform independent Bernoulli trials and count failures until $t$ successes occur.
  2. If $t$ is non-integer (including $t = 1/2$) with $p = 1/2$, generate a gamma deviate $X \sim \text{Gamma}(t, \theta)$ with $\theta = (1-p)/p$, then generate $N \sim \text{Poisson}(X)$.

Both methods yield integers $N$ with probabilities

$$ \Pr[N=n] = \binom{t+n-1}{n} p^t (1-p)^n, $$

as required.

∎

This solution corrects the critical error in the original solution by replacing the unjustified exponential–geometric mapping with the mathematically justified gamma–Poisson mixture approach.