When you cannot solve it, simulate it

The escape hatch

Analytic probability runs out quickly. The distribution of the maximum of three correlated variables, the chance that a queue exceeds capacity in the next hour, the spread of an A/B test result under a realistic model of user behaviour — all of these are awkward or impossible to write down in closed form.

All of them are easy to simulate. If you can describe how the randomness works as a procedure, you can run the procedure a hundred thousand times and count. This is called Monte Carlo, and it is the most practically useful idea in applied probability.

The pattern, in six lines

import numpy as np
rng = np.random.default_rng(0)

trials = 200_000
result = rng.normal(0, 1, trials) + rng.normal(0, 1, trials) > 1.5
print(result.mean())      # 0.1444

That estimates the probability that the sum of two standard normals exceeds 1.5. The exact answer is available here — the sum is normal with standard deviation √2, giving 0.1444 — which is exactly why it is a good first example: the simulation agrees, so you can trust the method on the problems where the exact answer is not available.

How accurate is the estimate

The estimate is a proportion from n independent trials, so its standard error is

se = sqrt( p (1 - p) / n )

With p ≈ 0.14 and n = 200,000:

se = sqrt(0.14 * 0.86 / 200000) = 0.00078

So the answer is 0.144 give or take about 0.0008 — three reliable digits. To get one more digit you need a hundred times more trials, because the error falls with √n. That relationship is the whole cost model of simulation: cheap to be roughly right, expensive to be very precise.

A problem worth simulating: the queue

Requests arrive at a service at 100 per second on average, and each takes 8 milliseconds. Average utilisation is 80 per cent, which sounds comfortable. How often does a request wait more than 100 ms?

The analytic answer needs queueing theory. The simulation needs a loop.

import numpy as np
rng = np.random.default_rng(1)
n = 200_000
gaps = rng.exponential(1/100, n)          # arrivals, seconds
service = rng.exponential(0.008, n)       # service times, seconds

arrival = np.cumsum(gaps)
finish = np.zeros(n)
prev = 0.0
for i in range(n):                        # a queue is inherently sequential
    start = max(arrival[i], prev)
    prev = finish[i] = start + service[i]
wait = finish - arrival - service
print((wait > 0.1).mean(), np.percentile(wait, 99))

You will find several per cent of requests waiting over 100 ms at 80 per cent utilisation, and a 99th percentile far above the 8 ms average. That gap between the mean and the tail is the single most important fact about capacity planning, and it emerges from twelve lines you can run on a laptop.

Two things that ruin a simulation

Correlations you did not model. If arrivals are bursty rather than independent — and real traffic always is — the simulation above understates the tail substantially. A simulation is only as good as the generating story you gave it, and a wrong story produces confident wrong numbers with no warning.

A bad random source, or a shared one. Use np.random.default_rng(seed) rather than the legacy global np.random.seed. The modern generator has better statistical properties, and passing an explicit generator object avoids the situation where two parts of your code silently share and reorder each other's random stream. When simulations must be reproducible, seed explicitly and record the seed with the results.

Sampling from distributions you have

Every library gives you the standard ones directly: rng.normal, rng.binomial, rng.poisson, rng.exponential, rng.choice. For a distribution you have only as a set of weights, rng.choice(values, p=weights) handles it, and that is exactly what a language model's sampler does at each step — draw one token from a categorical distribution over 50,000 options.

For resampling from data you already have, rather than a distribution you assumed, the operation is rng.choice(data, size=len(data), replace=True). That is the bootstrap, and it gets a full lesson in the statistics module.

Where simulation beats a formula even when a formula exists

Three cases:

• Communication. "In 200,000 simulated weeks, we breached capacity in 3 per cent of them" is understood by everyone in the room. A closed-form tail probability is not.
• Composition. Formulas do not compose; simulations do. Add a retry policy, a cache, a second server, and the formula must be re-derived while the simulation gains four lines.
• Checking your algebra. If you derived a result analytically, simulate it. Ten minutes, and it catches the sign error that would otherwise ship.

The rule to keep

If you can write down the procedure, you can estimate the probability. Accuracy improves with the square root of the trial count, so budget accordingly, and remember that the simulation tests your assumptions rather than the world.