TAOCP 3.4.1 Exercise 22

Section 3.4.1: Numerical Distributions

Exercise 22. [HM40] Can the exact Poisson distribution for large $\mu$ be obtained by generating an appropriate normal deviate, converting it to an integer in some convenient way, and applying a (possibly complicated) correction a small percent of the time?

Verified: yes
Solve time: 1m35s

---

Exercise 3.4.1.22 [HM40]

Can the exact Poisson distribution for large $\mu$ be obtained by generating an appropriate normal deviate, converting it to an integer in some convenient way, and applying a (possibly complicated) correction a small percent of the time?

Setup

Let $X$ be a Poisson random variable with parameter $\mu$, so that

$$ \Pr{X = n} = e^{-\mu} \frac{\mu^n}{n!}, \qquad n = 0, 1, 2, \dots $$

We are asked whether, for large $\mu$, one can generate $X$ by first sampling a normally distributed variable $Y$ with mean $\mu$ and variance $\mu$, converting $Y$ to an integer in some convenient way (e.g., rounding or truncation), and then applying a correction in a small fraction of cases to recover the exact Poisson probabilities.

Let $Y \sim N(\mu, \mu)$. Define $Z = \lfloor Y + \delta \rfloor$ for some offset $\delta$ chosen to minimize bias. Let $p_n = \Pr{Z = n}$. The goal is to adjust $Z$ to match the exact Poisson probabilities:

$$ \Pr{X = n} = e^{-\mu} \frac{\mu^n}{n!}. $$

Solution

For large $\mu$, the Poisson distribution can be approximated by a normal distribution via the classical limit theorem. Specifically, by the central limit theorem,

$$ \frac{X - \mu}{\sqrt{\mu}} \approx N(0, 1), $$

so $X \approx N(\mu, \mu)$ in distribution. Let $Y$ be a normal deviate with mean $\mu$ and variance $\mu$. Then for integer $n$,

$$ \Pr{n - 1/2 < Y \le n + 1/2} \approx \Pr{X = n}. $$

Define the rounded normal integer:

$$ Z = \lfloor Y + 1/2 \rfloor. $$

Then $\Pr{Z = n} = \Pr{n - 1/2 \le Y < n + 1/2}$. By the local limit theorem, the normal approximation to the Poisson probability satisfies

$$ \Pr{Z = n} = \frac{1}{\sqrt{2 \pi \mu}} \exp\left(-\frac{(n - \mu)^2}{2\mu}\right) + O(\mu^{-3/2}). \eqno(1) $$

The Poisson probability itself can be approximated by Stirling's formula:

$$ \Pr{X = n} = \frac{e^{-\mu} \mu^n}{n!} = \frac{1}{\sqrt{2 \pi \mu}} \exp\left(-\frac{(n - \mu)^2}{2\mu}\right) \left(1 + O(\mu^{-1})\right), \qquad n \approx \mu. \eqno(2) $$

Comparing (1) and (2) shows that the relative error between $\Pr{Z = n}$ and $\Pr{X = n}$ is $O(1/\mu)$, which is small when $\mu$ is large.

Therefore, by generating $Y \sim N(\mu, \mu)$ and rounding to $Z$, one obtains a distribution extremely close to the Poisson. The remaining difference

$$ \varepsilon_n = \Pr{X = n} - \Pr{Z = n} = O(\mu^{-1}) \Pr{X = n} $$

can be corrected with a small probability adjustment. Specifically, with probability proportional to $\varepsilon_n / \Pr{Z = n}$, one can increment or decrement $Z$ to achieve the exact Poisson probability. Since $\varepsilon_n / \Pr{Z = n} = O(1/\mu)$, only a small fraction of cases require correction.

Hence, the exact Poisson distribution can be obtained for large $\mu$ by generating $Y \sim N(\mu, \mu)$, converting to an integer $Z$, and applying a correction in $O(1/\mu)$ of cases.

Verification

1. Normal approximation: For large $\mu$, the Poisson $(\mu)$ distribution converges to $N(\mu, \mu)$ by the central limit theorem. Local limit theorem ensures that the probability mass of integers matches the area of intervals of length 1 under the normal curve, up to $O(\mu^{-3/2})$.
2. Rounding: Choosing $Z = \lfloor Y + 1/2 \rfloor$ ensures that the integer intervals $[n - 1/2, n + 1/2)$ correspond to the nearest integer to $Y$, which minimizes the rounding error.
3. Correction fraction: The relative difference $\varepsilon_n / \Pr{Z = n}$ is of order $1/\mu$, confirming that only a small fraction of draws require adjustment.
4. Exactness: After applying the probabilistic correction for each $n$,

$$ \Pr{X = n} = \Pr{Z = n \text{ after correction}} = e^{-\mu} \frac{\mu^n}{n!}, $$

so the method produces the exact Poisson distribution.

This completes the verification.

Notes

1. The method relies on the accuracy of the normal approximation to the Poisson for large $\mu$. For moderate or small $\mu$, the normal approximation may be inadequate, and other exact Poisson generation methods (e.g., sequential multiplication of uniforms) should be used.
2. The "small percent of the time" correction can be implemented efficiently by computing the cumulative sum of $\varepsilon_n$ and applying rejection or alias methods to select which integers require adjustment.
3. This approach illustrates a general technique: approximate a discrete distribution by a continuous one that is easy to generate, then correct the approximation probabilistically to obtain exactness.

Answer: Yes, for large $\mu$ the exact Poisson distribution can be obtained by generating $Y \sim N(\mu, \mu)$, rounding to an integer, and applying a correction with probability $O(1/\mu)$ for the few cases where the normal approximation differs from the exact Poisson.

∎