Probability by simulation

Probability questions are easy to get wrong by reasoning about them and hard to get wrong by counting. This section takes that seriously: for each idea you simulate the situation tens of thousands of times in a few lines of Python and count what happened, then do the maths and check the two agree. The habit is worth more than any single result here, because it means you can settle an argument about a probability on your own instead of looking for someone to confirm it.

Count what happens#

Every simulation in this section has the same three steps, and they are worth naming because the third is the one people leave out. Generate the situation many times. Reduce each repetition to a true or false. Take the mean of those booleans, which in NumPy is the share that were true. That is the probability. Two dice summing to seven is the smallest example: there are 36 equally likely pairs and six of them sum to seven, so the answer is one sixth, and the simulation should land on 0.167.

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

rolls = rng.integers(1, 7, size=(100_000, 2))       # two dice, 100,000 times
print((rolls.sum(axis=1) == 7).mean())              # 0.16775, against 6/36 = 0.1667

Any probability is the long-run share of times something happens. If you can simulate the situation, you can estimate it by counting, and check an exact answer you worked out on paper.

How wrong can a simulation be?#

A simulated answer is itself a measurement, so it has its own error bar, and knowing it is what separates using this technique from waving at it. The estimate is a proportion from n independent trials, so the margin is the same one as for any proportion: 1.96 times the square root of p times one minus p over n. At 100,000 trials and a true answer near one sixth, that is about two parts in a thousand, which is why 0.16775 and 0.1667 agreeing to two decimals is a pass and agreeing to four would be luck.

import math

def mc_margin(p_hat, trials):
    return 1.96 * math.sqrt(p_hat * (1 - p_hat) / trials)

for trials in (1_000, 100_000, 10_000_000):
    print(trials, round(mc_margin(0.1667, trials), 5))
# 1000 0.0231
# 100000 0.00231
# 10000000 0.00023

More samples, closer answer#

That square-root law is the law of large numbers with a rate attached. The average of repeated draws converges on the expectation, but slowly and not monotonically: it can be further away after more draws than it was before. Watch it wander in on a fair coin.

for n in [10, 100, 1_000, 100_000]:
    print(n, rng.integers(0, 2, n).mean())   # coin flips: drifts towards 0.5
# 10 0.8
# 100 0.49
# 1000 0.509
# 100000 0.49946

The birthday problem#

Twenty-three people in a room, and the chance that two share a birthday is better than even. Nearly everyone guesses far lower, because they answer a different question — the chance that someone shares their birthday — when the real one counts all 253 pairs. It is the standard demonstration that intuition about compound events is not to be trusted, and it has an exact answer to check against: multiply 365/365 by 364/365 by 363/365 and so on, which is the probability that everyone is different, then subtract from one.

def shared_birthday(people, trials=50_000):
    days = rng.integers(0, 365, size=(trials, people))
    s = np.sort(days, axis=1)
    return (np.diff(s, axis=1) == 0).any(axis=1).mean()

def exact(people):
    none_shared = 1.0
    for i in range(people):
        none_shared *= (365 - i) / 365
    return 1 - none_shared

for k in (10, 23, 40, 60):
    print(k, round(shared_birthday(k), 3), round(exact(k), 3))
# 10 0.116 0.117
# 23 0.509 0.507
# 40 0.893 0.891
# 60 0.994 0.994

The engineering version of this is hash collisions, and the numbers are the same shape. With a 64-bit hash you expect a collision after about four billion keys rather than after eighteen quintillion, because it is pairs that matter, not items. Any time you reason about the chance that some pair of things collides, the birthday calculation is the one to reach for.

Monty Hall#

Three doors, a car behind one, you pick a door, the host — who knows where the car is — opens a different door showing a goat and offers you the switch. Almost everyone says it makes no difference. It doubles your chances, and the simulation is shorter than the argument. The key is that your first pick is right one time in three, and the host's behaviour cannot change that, so the other unopened door holds the car the remaining two times in three. Switching wins exactly when your first pick was wrong.

trials = 200_000
car = rng.integers(0, 3, trials)
first_pick = rng.integers(0, 3, trials)

print("stay  ", (first_pick == car).mean())     # 0.332165
print("switch", (first_pick != car).mean())     # 0.667835

Base rates#

A test catches 99% of people with a condition and wrongly flags 5% of people without it. The condition affects 1% of people. You test positive. How likely is it you have it?

n = 1_000_000
sick = rng.random(n) < 0.01
positive = np.where(sick, rng.random(n) < 0.99, rng.random(n) < 0.05)
print(sick[positive].mean())    # 0.16511298003850441: most positives are false alarms

Because healthy people vastly outnumber sick ones, 5% of them outnumber 99% of the sick. That is Bayes' theorem done in counts rather than symbols, and the same arithmetic explains why a fraud model with a low false positive rate still produces mostly false alerts.

Getting there by simulation rather than by formula is not a lesser route. The three lines above encode the story directly — who is sick, what the test does to each group — so a misreading of the question shows up as code you cannot write, whereas a misreading of the formula shows up as a confident wrong number. When a conditional probability argument matters, write both and make them agree.

Averages become bell curves#

Waiting times are about as un-bell-shaped as a quantity gets: most waits are short, a few are very long, and the histogram is a cliff falling to the right. Take thirty of them and average, though, and something strange happens to the average. Its distribution is close to symmetric, close to normal, and much narrower than the waits themselves.

waits = rng.exponential(scale=2.0, size=(100_000, 30))   # skewed, not bell-shaped
means = waits.mean(axis=1)                                 # the average of 30 waits
print(means.mean(), means.std())
# 2.00171649553826 0.3669219319476103
# the mean of the averages is the mean of the waits; 2.0 / sqrt(30) = 0.365

The central limit theorem is the claim that this happens for almost anything you average, whatever shape it started as, and that the limit is specifically the normal distribution. Claims like that deserve measurement rather than assent, so measure it. Standardise the average — subtract its mean, multiply by the square root of n — and a normal result must have its middle 95 per cent between minus 1.96 and plus 1.96, and zero skewness. Then watch those three numbers as n grows.

for n in (1, 2, 5, 30, 200):
    m = rng.exponential(1.0, size=(200_000, n)).mean(axis=1)
    z = (m - 1.0) * math.sqrt(n)          # standardised: normal means 2.5% at -1.96
    print(f"n={n:4d}  2.5% {np.percentile(z, 2.5):6.2f}   97.5% {np.percentile(z, 97.5):5.2f}"
          f"   skew {float(((z - z.mean()) ** 3).mean()):5.2f}")

# n=   1  2.5%  -0.97   97.5%  2.68   skew  2.00
# n=   2  2.5%  -1.24   97.5%  2.51   skew  1.36
# n=   5  2.5%  -1.51   97.5%  2.35   skew  0.90
# n=  30  2.5%  -1.78   97.5%  2.12   skew  0.37
# n= 200  2.5%  -1.89   97.5%  2.03   skew  0.14
# normal reference:  -1.96, 1.96, skew 0

That last point is the practical one, and it is why this is worth measuring rather than quoting. The average of many independent values is approximately normally distributed whatever the shape of the original values, with spread shrinking as the square root of the sample size, which is why error bars on averages work. But approximately is doing real work in that sentence. The rule of thumb that thirty is enough comes from mildly skewed data; on revenue per user, or latency, or anything where one observation in a thousand is a hundred times the others, thirty is nowhere near enough and the error bar you computed is too narrow. Simulate your own data before trusting a normal interval on it.