How to Generate a Poisson Random Variable from a Uniform One

A single standard uniform random variable can be transformed into a Poisson random variable by dividing the interval from 0 to 1 into pieces whose lengths match the Poisson probabilities. This is the inverse-transform method. It gives a direct algorithm for one sample; to generate independent Poisson samples, use independent uniform draws.

The inverse-transform idea

Let $U$ be uniform on $[0,1)$, and let $X$ have a Poisson distribution with rate parameter $\lambda\geq 0$. Its probability mass function is

$$P(X=k)=p_k=e^{-\lambda}\frac{\lambda^k}{k!},\qquad k=0,1,2,\ldots$$

Define the cumulative distribution function by $F(k)=P(X\leq k)=\sum_{j=0}^{k}p_j$. The inverse-transform rule is to return the smallest nonnegative integer $k$ for which $F(k)\geq U$.

Equivalently, assign the interval from $F(k-1)$ to $F(k)$ to the outcome $k$, with $F(-1)=0$. The interval’s length is $F(k)-F(k-1)=p_k$. A uniform draw lands in that interval with probability $p_k$, exactly the probability required for a Poisson outcome. The precise assignment of interval endpoints does not affect the distribution because any individual endpoint has probability zero.

A step-by-step algorithm

You do not need factorials or a table of probabilities. Start with $p_0=e^{-\lambda}$, then use the recurrence

$$p_k=p_{k-1}\frac{\lambda}{k},\qquad k\geq 1.$$

Maintain a running cumulative total. Stop as soon as that total is at least $U$.

1. Set $k=0$, $p=e^{-\lambda}$, and $C=p$.
2. While $U>C$, increase $k$ by 1, update $p\leftarrow p\lambda/k$, and add $p$ to $C$.
3. Return $k$.

For $\lambda=0$, return 0 directly. In an ideal mathematical calculation, the cumulative probabilities tend to 1 as $k$ increases. A computer-generated uniform is normally in $[0,1)$, so a finite count is returned.

Worked example

Take $\lambda=3$ and suppose the uniform draw is $U=0.8$. The probability at zero is $p_0=e^{-3}\approx 0.049787$. The recurrence gives the following probabilities and cumulative totals:

$k$	$p_k$ (approximately)	$F(k)$ (approximately)
0	0.049787	0.049787
1	0.149361	0.199148
2	0.224042	0.423190
3	0.224042	0.647232
4	0.168031	0.815263

The cumulative total is below 0.8 at $k=3$ and above it at $k=4$. Therefore the transformed value is $X=4$. In general, values of $U$ between about 0.647232 and 0.815263 produce 4 for this rate.

Python implementation

The function below transforms a supplied uniform value. It uses only the Python standard library’s math module and is suitable for ordinary, moderate values of $\lambda$.

import math

def poisson_from_uniform(u, lam):
    if not 0.0 <= u < 1.0:
        raise ValueError('u must be in [0, 1)')
    if lam < 0.0:
        raise ValueError('lam must be nonnegative')
    if lam == 0.0:
        return 0

    k = 0
    p = math.exp(-lam)
    cumulative = p

    while u > cumulative:
        k += 1
        p *= lam / k
        cumulative += p

    return k

For example, poisson_from_uniform(0.8, 3.0) returns 4. To generate a sample in a program, pass in a uniform value from a random-number generator; repeated samples intended to be independent should use independent draws.

Why the algorithm is correct

For a given $k$, the procedure returns $k$ when the uniform draw falls in the portion of the unit interval allocated to that count. That portion has length $F(k)-F(k-1)=p_k$. Since $U$ is uniform, the probability of landing there is its length, so $P(X=k)=p_k$. This verifies the full distribution, not just its mean or variance.

Common mistakes and numerical limits

• Stopping at the first probability above $U$ instead of the cumulative total. Compare $U$ with $F(k)$, the sum through $k$, not with the individual probability $p_k$.
• Using the wrong recurrence index. When advancing from $k-1$ to $k$, multiply by $\lambda/k$.
• Treating $\lambda$ as an integer. The rate can be any nonnegative real number; only the generated outcome must be an integer.
• Passing $U=1$. For $\lambda>0$, the inverse CDF does not reach 1 at any finite count in exact arithmetic; for $\lambda=0$, it reaches 1 at $k=0$. Standard pseudorandom generators ordinarily return values below 1; if your source can return 1, handle that endpoint separately.
• Using this direct recurrence for very large $\lambda$. It can require many iterations, and computing $e^{-\lambda}$ may underflow to zero in floating-point arithmetic. For large rates, use a tested Poisson random-variate routine based on a method designed for that range rather than relying on this simple inverse-CDF implementation.

The inverse-transform construction itself applies to any discrete distribution: use its cumulative probabilities and return the first outcome whose cumulative total reaches the uniform draw.