OS-015 · Stochastic processes

Poisson Process Simulation

A Poisson process scatters events in time with no memory and no pattern: calls to a help line, clicks of a Geiger counter, arrivals at a queue. Its counts follow the Poisson distribution, its gaps are exponential, and the gap you happen to land in is, on average, twice as long as a typical one.

Interactive model

OS-015 · poisson-process · seeded
eventsgap containing Thorizon T
runs with k eventsPoisson(λT) expectation
gaps (density)λe^(−λx)mean 1/λ
Parameters
Independent realizations of the process on [0, T]
Same seed, same run — here, in simulation.js and in simulation.py.
Mean count per run
20.053
Variance of counts
18.220
Theory for both: λT
20.000
Mean gap
0.4984
Theory: 1/λ
0.5000
Gap SD / mean (theory 1)
1.003
Mean gap containing T
1.0085
Theory: (2 − e^(−λT))/λ
1.0000

Example run: default parameters, seed 42.

What the simulation shows

A homogeneous Poisson process with rate λ has three equivalent descriptions, and the lab lets you check each against the others:

  • Gaps. The times between consecutive events are independent and exponentially distributed with mean 1/λ. This is how the simulation generates the process.
  • Counts. The number of events in any interval of length T has the Poisson distribution with mean λT, and counts in non-overlapping intervals are independent.
  • Small intervals. In a short interval of length h an event occurs with probability λh + o(h), two or more with probability o(h), independently of the past.

Each run generates events on [0, T]. The timeline shows the latest runs; the histograms collect all of them. Press Step to add one run at a time and watch the count histogram fill in toward the Poisson expectation drawn over it.

  1. Set t = 0.
  2. Draw u ~ U[0, 1) and set the gap g = −ln(1 − u)/λ (inversion of the exponential distribution function).
  3. Add g to t. If t > T, stop; otherwise record an event at t and go to step 2.

Why the counts are Poisson

Split [0, T] into n tiny slots, each holding an event with probability λT/n. The count is binomial, and as n → ∞ the binomial converges to

P(N(T) = k) = e−λT (λT)k / k!

This “law of rare events” is why the Poisson distribution appears whenever many independent opportunities each have a small chance. Its mean and variance are both λT, so their ratio, the index of dispersion, is 1. The lab reports the sample mean and variance of the counts next to λT; real count data with variance well above the mean are a sign of clustering or a varying rate, which the model excludes.

A second way to simulate the process uses a further property: given that N(T) = n, the n event times are distributed like n independent uniform points on [0, T], sorted. So one can draw n from the Poisson distribution first and then scatter uniform points. Both methods produce the same process.

Memoryless gaps and the inspection paradox

The exponential distribution is memoryless: having waited s for the next event, the remaining wait is again exponential with the same mean, P(G > s + t | G > s) = P(G > t). It is the only continuous distribution with this property, and it is what makes the process “forget” its past.

Now look at the gap that contains the fixed time T, shaded in the timeline. Long gaps cover more of the time axis, so T is more likely to land in one of them. By memorylessness the wait from T to the next event is exponential with mean 1/λ, and the time back to the previous event (or to 0) has mean (1 − e−λT)/λ, so

E[gap containing T] = (2 − e−λT) / λ ≈ 2/λ

This is the waiting-time paradox: buses arriving as a Poisson process every 10 minutes on average leave a randomly arriving passenger in a gap that averages 20 minutes, and the passenger still waits 10 on average.

The same effect decides how the gaps must be averaged. The lab averages every gap it draws, including the one that crosses T. Keeping only gaps that end before T would favor short gaps and bias the mean below 1/λ, by roughly 1/(λ²T).

Example result

With λ = 2 events per unit time, T = 10 and 1,000 runs (seed 42), the counts have mean 20.053 and variance 18.220, both estimates of λT = 20.000; the variance is the noisier estimate, with a standard error of about λT·√(2/R) ≈ 0.9 here. The gaps average 0.4984 against 1/λ = 0.5000, with a ratio of standard deviation to mean of 1.003, as an exponential distribution requires. The gap containing T averaged 1.0085, close to the predicted 1.0000: twice a typical gap.

Output of poisson-process with rate=2 horizon=10 runs=1000 seed=42
Mean count per run20.053
Variance of counts18.220
Theory for both: λT20.000
Mean gap0.4984
Theory: 1/λ0.5000
Gap SD / mean (theory 1)1.003
Mean gap containing T1.0085
Theory: (2 − e^(−λT))/λ1.0000

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

JavaScript — simulation.js
// One realization on [0, T]: add exponential gaps until the clock passes T.
// Returns the event times and every gap drawn, including the last one, which
// straddles T (its event falls after the horizon and is not counted).
export function simulateRun(rate, T, rng) {
  const times = [], gaps = [];
  let t = 0;
  for (;;) {
    const g = rng.exponential(rate); // inter-arrival time ~ Exp(λ), mean 1/λ
    gaps.push(g);
    t += g;
    if (t > T) return { times, gaps };
    times.push(t);
  }
}

// Running totals over many runs: counts N(T) and inter-arrival gaps.
export function createStats() {
  return { runs: 0, sumN: 0, sumN2: 0, gapN: 0, sumG: 0, sumG2: 0, sumStraddle: 0, hist: [] };
}
export function addRun(s, r) {
  const n = r.times.length;
  s.runs++; s.sumN += n; s.sumN2 += n * n;
  s.hist[n] = (s.hist[n] || 0) + 1;
  for (const g of r.gaps) { s.gapN++; s.sumG += g; s.sumG2 += g * g; }
  s.sumStraddle += r.gaps[r.gaps.length - 1];
}
Python — simulation.py
def simulate_run(rate, T, rng):
    """Event times on [0, T] from exponential gaps; the last gap drawn straddles T."""
    times, gaps = [], []
    t = 0.0
    while True:
        g = rng.exponential(rate)  # inter-arrival time ~ Exp(rate), mean 1/rate
        gaps.append(g)
        t += g
        if t > T:
            return times, gaps
        times.append(t)

Run the full model from a terminal:

Shell
node simulation.js --rate 2 --horizon 10 --seed 42
python simulation.py --rate 2 --horizon 10 --seed 42

Download the model

poisson-process.zip6 files · 9.2 KBDownload .zip
Files in poisson-process.zip
FileContentsSize
README.mdHow to run it, parameters, example result2.2 KB
model.jsonParameters with ranges and defaults, outputs, the seed-42 example2.0 KB
simulation.jsThe model used on this page plus a command-line runner (Node 18+)7.6 KB
simulation.pyPython 3.9+ port, standard library only6.2 KB
example.csvResult table from the example run899 B
LICENSE.txtMIT License1.0 KB

Assumptions & limitations

  • Constant rate. The rate does not change with time of day or season. A time-varying rate gives an inhomogeneous Poisson process, usually simulated by thinning a faster homogeneous one.
  • Independence. An event neither triggers nor suppresses others. Aftershocks, retweets and bursts of traffic violate this; self-exciting (Hawkes) processes model them.
  • No simultaneous events. Two events never occur at exactly the same time (in exact arithmetic).
  • Arrivals only. The process describes when things happen, not how long they take; queues add service times on top, as in the M/M/1 queue.

References

  1. Kingman, J. F. C. (1993). Poisson Processes. Oxford University Press.
  2. Feller, W. (1968). An Introduction to Probability Theory and Its Applications, Vol. 1, 3rd ed. Wiley.
  3. Grinstead, C. M. & Snell, J. L. Introduction to Probability, chapter 5 (important distributions, including the Poisson). American Mathematical Society. math.dartmouth.edu/~prob/prob/prob.pdf