51. Poisson processes

🎬 Video · 15 min
💡 Every code box below is live — edit it and hit Run.

A Markov chain moves in discrete steps. A Poisson process models events arriving in continuous time — customers, phone calls, radioactive decays, website hits.

It is the canonical model for "things happen at random, at a steady average rate", and it ties together three distributions you already know.

Three equivalent definitions

A counting process N(t) — the number of events by time t — is a Poisson process with rate \lambda if any of the following hold. They're equivalent, and each is useful for different purposes.

1. By the counts. N(0) = 0, increments over disjoint intervals are independent, and

N(t) \sim \text{Poisson}(\lambda t)

2. By the gaps. The waiting times between consecutive events are independent \text{Exponential}(\lambda) variables.

3. Infinitesimally. In a short interval of length h:

P(\text{1 event}) = \lambda h + o(h), \quad P(\ge 2 \text{ events}) = o(h)

with disjoint intervals independent.

The second definition is the one to hold in your head, because it makes simulation trivial and explains the memorylessness. The first is what you compute with.

Why they're the same

The link between them is a one-line calculation.

"No events by time t" means the first waiting time exceeds t:

P(N(t) = 0) = P(T_1 > t) = e^{-\lambda t}

And that's exactly the Poisson probability of zero events with mean \lambda t:

P(N(t) = 0) = \frac{(\lambda t)^0 e^{-\lambda t}}{0!} = e^{-\lambda t} \quad\checkmark

So Exponential gaps and Poisson counts are two descriptions of one process. The Exponential distribution's memorylessness (§1) is what makes the third definition work: the chance of an event in the next instant never depends on how long you've waited.

Key properties

Mean and variance: E[N(t)] = \operatorname{Var}(N(t)) = \lambda t. Equal, as for any Poisson.

Superposition. Merge two independent Poisson processes with rates \lambda_1 and \lambda_2, and you get a Poisson process with rate \lambda_1 + \lambda_2. Rates add.

Thinning. Keep each event independently with probability p, and the result is Poisson with rate \lambda p — and the kept and discarded streams are independent Poisson processes, which is genuinely surprising.

Conditional uniformity. Given that exactly n events occurred in [0, T], their times are distributed as n independent uniform points on [0, T] (order statistics, §2). The process has no preferred moments.

That last property is the precise sense in which a Poisson process is "completely random", and it's how you simulate a fixed number of events.

The clustering illusion

Uniform randomness produces clumps, and people consistently misread them.

Scatter 100 points uniformly on a line and you'll see visible clusters and visible gaps. That's not evidence of an underlying cause — it's what randomness looks like. Genuinely evenly spaced points would be far less likely.

This is why "cancer clusters", bombing patterns, and apparent runs of bad luck so often survive scrutiny as coincidence. The test isn't "does it look clumped?" but "is it more clumped than a Poisson process would be?"

When it doesn't apply

The assumptions are strong:

  • Constant rate. Real arrivals usually vary by hour and day. Fix: an inhomogeneous Poisson process with \lambda(t).
  • Independence. Events that trigger further events — earthquake aftershocks, viral shares — cluster far more than Poisson allows. Fix: Hawkes or Cox processes.
  • No simultaneity. If arrivals come in batches (a bus full of customers), use a compound Poisson process.

The diagnostic is the same as §6's overdispersion check: if the variance of counts exceeds their mean, the process isn't Poisson.

Worked example

A call centre receives calls at \lambda = 5 per hour.

(a) Probability of exactly 3 calls in one hour.

N(1) \sim \text{Poisson}(5):

P(N = 3) = \frac{5^3e^{-5}}{3!} = \frac{125 \times 0.006738}{6} \approx 0.1404

(b) Probability of no calls in 30 minutes.

Now t = 0.5, so the mean is \lambda t = 2.5:

P(N(0.5) = 0) = e^{-2.5} \approx 0.0821

(c) Expected time until the next call.

Gaps are \text{Exponential}(5):

E[T] = \frac{1}{5} = 0.2 \text{ hours} = 12 \text{ minutes}

(d) You've already waited 20 minutes. Expected additional wait?

Still 12 minutes — the Exponential is memoryless. The process has no obligation to "catch up".

(e) Given 10 calls arrived between 9am and 10am, when did they arrive?

By conditional uniformity, the 10 arrival times are distributed as 10 independent uniform draws on that hour. Knowing the count tells you nothing about the pattern.

Doing it in Python

The two equivalent constructions:

import numpy as np
from scipy.stats import poisson, expon

rng = np.random.default_rng(0)
lam, T, trials = 5.0, 1.0, 200_000

# Construction 1: draw the count directly
counts_direct = rng.poisson(lam * T, trials)

# Construction 2: accumulate exponential gaps until we pass T
def count_via_gaps(lam, T, rng):
    t, n = 0.0, 0
    while True:
        t += rng.exponential(1 / lam)
        if t > T:
            return n
        n += 1

counts_gaps = np.array([count_via_gaps(lam, T, rng) for _ in range(20_000)])

print(f"direct  : mean {counts_direct.mean():.4f}  var {counts_direct.var():.4f}")
print(f"via gaps: mean {counts_gaps.mean():.4f}  var {counts_gaps.var():.4f}")
print(f"theory  : mean {lam*T}       var {lam*T}")
print("\nSame process, two descriptions.")

The call centre questions:

import numpy as np
from scipy.stats import poisson

lam = 5.0

print(f"(a) P(3 calls in 1 hour)   = {poisson.pmf(3, lam * 1.0):.4f}")
print(f"(b) P(0 calls in 30 min)   = {poisson.pmf(0, lam * 0.5):.4f}")
print(f"    = e^-2.5               = {np.exp(-2.5):.4f}")
print(f"(c) expected wait          = {1/lam:.4f} hours = {60/lam:.1f} minutes")
print(f"(d) after waiting 20 min   = {60/lam:.1f} minutes (memoryless)")

Superposition and thinning:

import numpy as np

rng = np.random.default_rng(1)
T, trials = 1.0, 200_000

a = rng.poisson(3.0 * T, trials)          # rate 3
b = rng.poisson(2.0 * T, trials)          # rate 2

print("superposition: merging rate-3 and rate-2 streams")
print(f"  merged mean {(a+b).mean():.4f}  var {(a+b).var():.4f}   (theory 5, 5)")

# Thinning: keep each event with probability 0.3
total = rng.poisson(10.0 * T, trials)
kept = rng.binomial(total, 0.3)
discarded = total - kept

print("\nthinning a rate-10 stream with p = 0.3")
print(f"  kept      mean {kept.mean():.4f}  var {kept.var():.4f}   (theory 3, 3)")
print(f"  discarded mean {discarded.mean():.4f}  var {discarded.var():.4f}   (theory 7, 7)")
print(f"  correlation between kept and discarded: "
      f"{np.corrcoef(kept, discarded)[0,1]:+.4f}   <- independent!")

Conditional uniformity, and the clustering illusion:

import numpy as np

rng = np.random.default_rng(2)

# Given n events in [0, T], their times are uniform
n, T = 12, 1.0
times = np.sort(rng.uniform(0, T, n))
print("12 arrival times in one hour:", (times * 60).round(1))

gaps = np.diff(np.r_[0, times, T]) * 60
print(f"gaps (minutes): {gaps.round(1)}")
print(f"  shortest {gaps.min():.1f}, longest {gaps.max():.1f} "
      f"-- evenly spaced would be {60/(n+1):.1f} each")

# How often does 'random' produce a visible cluster?
close_pairs = 0
for _ in range(20_000):
    t = np.sort(rng.uniform(0, 60, 12))
    close_pairs += (np.diff(t) < 1.0).any()       # two arrivals within a minute
print(f"\nP(some two of 12 arrivals fall within 1 minute) = {close_pairs/20_000:.4f}")
print("Clumping is the norm, not a signal.")

Detecting a process that isn't Poisson:

import numpy as np

rng = np.random.default_rng(3)
trials = 20_000

# Genuine Poisson
pois = rng.poisson(5.0, trials)

# Clustered arrivals: each 'trigger' brings a small batch
triggers = rng.poisson(1.5, trials)
clustered = np.array([rng.poisson(3.3, t).sum() if t else 0 for t in triggers])

for name, data in [("Poisson", pois), ("clustered", clustered)]:
    print(f"{name:>12}: mean {data.mean():.3f}  var {data.var():.3f}  "
          f"var/mean {data.var()/data.mean():.3f}")

print("\nvar/mean near 1 means Poisson. Well above 1 means overdispersed --")
print("the arrivals are clustering more than pure randomness allows.")

Your turn

1. Emails arrive at 12/hour. What's the probability of exactly 2 in 10 minutes?

2. Two independent Poisson processes have rates 3 and 7. What's the distribution of the merged process over one unit of time?

3. A machine fails as a Poisson process with rate 0.1/day. It's been running 30 days without failure. What's the expected additional time to failure?

Solutions

1. Ten minutes is t = 1/6 hour, so the mean is \lambda t = 12/6 = 2:

P(N = 2) = \frac{2^2e^{-2}}{2!} = \frac{4 \times 0.1353}{2} \approx 0.2707

About 27%. Note the whole calculation is "convert the window to a mean, then use Poisson" — the rate and the time only ever appear as their product.

2. By superposition, rates add:

N_1(t) + N_2(t) \sim \text{Poisson}\big((3 + 7)t\big)

Over one unit of time that's \text{Poisson}(10), with mean and variance both 10.

This also follows from §1's MGF argument: the sum of independent Poissons is Poisson with the summed rates. Superposition is that fact, restated for processes.

3. 10 days — the same as it was on day 1.

E[T] = \frac{1}{\lambda} = \frac{1}{0.1} = 10 \text{ days}

The 30 failure-free days are irrelevant, by memorylessness:

P(T > 30 + s \mid T > 30) = P(T > s)

This is the point at which the model should be questioned rather than trusted. Real machines do age — wear accumulates, and a 30-day-old machine is usually more likely to fail soon, not equally likely. A Poisson process assumes no ageing at all.

Where the assumption genuinely holds is in processes with no memory in the physics: radioactive decay, photon arrivals, and (approximately) calls arriving from a large independent population. For components that wear out, a Weibull lifetime with an increasing hazard rate is the standard alternative.

Since the course names it twice without defining it: the Weibull has survival function P(T > t) = e^{-(t/\lambda)^k} and hazard rate h(t) = \frac{k}{\lambda}\left(\frac{t}{\lambda}\right)^{k-1}. The shape parameter k is the whole point — k = 1 gives a constant hazard and collapses back to the Exponential exactly, k > 1 gives a hazard that rises with age (wear-out), and k < 1 a falling one (infant mortality, where surviving early failures makes a component more reliable). It is the Exponential with the memorylessness assumption relaxed in a single parameter, which is why it is the default in reliability engineering.

Check yourself in code

Verify that the two constructions of a Poisson process agree, and confirm that merging two processes adds their rates.

Print exactly this:

direct mean 5.0
gaps mean 5.0
merged mean 5.0
merged var 5.0
all match theory: True

Use default_rng(0) with 200000 direct draws at rate 5, 20000 gap-based simulations, and 200000 merged draws from rates 3 and 2. Round every value to 2 decimal places, and report True if all three means are within 0.1 of 5.

import numpy as np

rng = np.random.default_rng(0)
lam, T = 5.0, 1.0

direct = rng.poisson(lam * T, 200_000)
print("direct mean", round(direct.mean(), 2))

def count_via_gaps(lam, T, rng):
    t, n = 0.0, 0
    while True:
        t += rng.exponential(1 / lam)
        if t > T:
            return n
        n += 1

# Simulate 20000 counts via exponential gaps, then merge rate-3 and rate-2
# streams (200000 draws each) and report the merged mean and variance.
import numpy as np

rng = np.random.default_rng(0)
lam, T = 5.0, 1.0

direct = rng.poisson(lam * T, 200_000)
print("direct mean", round(direct.mean(), 2))

def count_via_gaps(lam, T, rng):
    t, n = 0.0, 0
    while True:
        t += rng.exponential(1 / lam)
        if t > T:
            return n
        n += 1

gaps = np.array([count_via_gaps(lam, T, rng) for _ in range(20_000)])
print("gaps mean", round(gaps.mean(), 2))

merged = rng.poisson(3.0 * T, 200_000) + rng.poisson(2.0 * T, 200_000)
print("merged mean", round(merged.mean(), 2))
print("merged var", round(merged.var(), 2))

means = [direct.mean(), gaps.mean(), merged.mean()]
print("all match theory:", all(abs(m - 5.0) < 0.1 for m in means))

A Poisson process is the continuous-time model of events at a steady random rate. Counts over an interval are Poisson(\lambda t); the gaps between events are independent Exponentials; and given the count, the arrival times are uniform. Rates add under superposition and scale under thinning. The assumptions — constant rate, independence, no batches — are strong, and the variance-to-mean ratio is how you check them.

Next: a process that moves rather than counts.