Poisson Distribution and Poisson Process in Python

Data · Published 2026-10-27

The Poisson distribution counts how many events happen in a fixed window when events arrive independently at a constant average rate λ. In Python, scipy.stats.poisson(mu=λ) gives you the probabilities, and a Poisson process (the arrival times themselves) is simulated by adding up exponential gaps with mean 1/λ.

Typical examples are support tickets per hour, requests per second on a server, typos per page or defects per batch. The examples were run with SciPy 1.18 and NumPy 2.5.

The distribution in one paragraph

If X ~ Poisson(λ), the probability of exactly k events is P(X = k) = λ^k e^−λ / k!, for k = 0, 1, 2, … Its single parameter λ is both the mean and the variance. That equality is the quickest diagnostic: if your counts have a variance much larger than their mean (overdispersion), a plain Poisson model will understate uncertainty, and a negative binomial is usually a better fit.

Probabilities with scipy.stats.poisson

Suppose a help desk receives on average 4 tickets per hour.

from scipy import stats

lam = 4                      # average events per hour
X = stats.poisson(mu=lam)
print(X.pmf([0, 2, 4, 6]).round(4))
print(round(X.cdf(2), 4))       # P(X <= 2)
print(round(X.sf(6), 4))        # P(X > 6)
print(X.mean(), X.var())
print(X.ppf(0.95))

[0.0183 0.1465 0.1954 0.1042]
0.2381
0.1107
4.0 4.0
8.0

Reading the output:

A completely quiet hour has a probability of about 1.8% (e⁻⁴).
Exactly four tickets, the most likely count together with three, happens about 19.5% of the time.
cdf(2) is the chance of two tickets or fewer: about 24%.
sf(6), the survival function, is P(X > 6) = 1 − cdf(6). Use sf rather than 1 - cdf: it stays accurate for tiny tail probabilities. Watch the boundary: P(X ≥ 7) is sf(6).
ppf(0.95) returns 8: in 95% of hours you will see at most 8 tickets, a natural staffing threshold.

Scaling the window

Rates scale with the window length. If the rate is 4 per hour, the count in a 15-minute window is Poisson(1) and the count over an 8-hour shift is Poisson(32). Always convert λ to the same unit as the window you ask about.

From counts to arrival times: the Poisson process

A homogeneous Poisson process with rate λ is the random sequence of event times behind those counts. It has three equivalent descriptions:

The number of events in any window of length t is Poisson(λt).
Counts in non-overlapping windows are independent.
The waiting times between consecutive events are independent exponential variables with mean 1/λ.

The third description gives the simplest simulation: draw exponential gaps and take their cumulative sum.

import numpy as np

rng = np.random.default_rng(7)
rate, horizon = 4.0, 1_000.0            # 4 events/hour over 1,000 hours
gaps = rng.exponential(scale=1 / rate, size=int(rate * horizon * 1.2))
times = np.cumsum(gaps)
times = times[times < horizon]
print(len(times))
print(times[:5].round(3))

4009
[0.177 0.433 0.575 0.799 0.851]

NumPy's exponential takes the scale (the mean gap, 1/λ), not the rate. Drawing 20% more gaps than the expected number makes it practically certain that the cumulative sum passes the horizon; the extra events are then cut off. Over 1,000 hours we expect 4,000 events and got 4,009.

Checking the simulation against the theory

Bin the arrival times into one-hour windows and compare the counts with the Poisson probabilities.

counts = np.histogram(times, bins=np.arange(0, horizon + 1, 1.0))[0]
print(round(counts.mean(), 3), round(counts.var(), 3))

obs = np.bincount(counts, minlength=11)[:11]
exp = X.pmf(np.arange(11)) * len(counts)
for k in range(0, 9, 2):
    print(k, obs[k], round(exp[k], 1))

4.009 3.931
0 18 18.3
2 144 146.5
4 187 195.4
6 110 104.2
8 33 29.8

The hourly mean and variance are both close to 4, and the observed frequencies track the expected ones. This is the counting description emerging from the waiting-time description.

An alternative simulation: count first, then place

Conditioned on the total count, the event times of a homogeneous Poisson process are uniformly scattered over the window. So you can draw the total from a Poisson distribution and then place that many points uniformly:

n = rng.poisson(rate * horizon)
times2 = np.sort(rng.uniform(0, horizon, size=n))
print(n)
print(round(np.diff(times2).mean(), 4))

3950
0.2532

The average gap is close to 1/λ = 0.25 hours, as it should be. This version is vectorised in one step and is convenient when you need a fixed window. The link with the continuous uniform distribution is not a coincidence: it is a defining property of the process.

Waiting-time questions

Because the gaps are exponential, waiting-time questions use scipy.stats.expon with scale=1/λ. The chance that no ticket arrives in the next 30 minutes:

print(round(stats.expon(scale=1 / rate).sf(0.5), 4))

0.1353

That is e⁻², identical to stats.poisson(2).pmf(0): "no event in half an hour" and "zero events in a Poisson(2) window" are the same statement. The exponential is also memoryless: having already waited 20 minutes does not change the distribution of the remaining wait.

When the Poisson model breaks down

Rates that vary over time. Traffic is busier at noon than at 3 a.m. Model this as a non-homogeneous process, for example by thinning: simulate at the peak rate and keep each event with probability λ(t) / λ_max. Counting by hour with pandas groupby is a quick way to see whether the rate is stable.
Clustering. Retries, cascading failures and viral bursts make events depend on each other. Variance then exceeds the mean.
A hard cap. If there are only n opportunities, each with probability p, the count is binomial. Poisson(np) is a good approximation only when n is large and p small.
Too many zeros. Data from many inactive units plus a few active ones needs a zero-inflated model.

Summary

Use stats.poisson(mu=λ) for count probabilities: pmf for exact counts, cdf and sf for tails, ppf for thresholds. Simulate a Poisson process with cumulative sums of rng.exponential(scale=1/λ), or by drawing a Poisson total and placing points uniformly. Check that the mean and variance of observed counts are similar before trusting the model.