Monte Carlo simulation

Why this matters

Many real probability questions have no tidy formula. What is the chance a project with 12 uncertain tasks finishes on time? How likely is a portfolio to lose more than 10% in a year? What share of manufactured assemblies will be out of tolerance when five parts each vary a little? Monte Carlo simulation answers all of them the same way: build a model of the random process, run it thousands of times, and count. It is a daily tool for data scientists, actuaries and operations research analysts, and it is also the best way to check your own formulas.

The recipe

1. Model one run of the process with random draws.
2. Repeat it N times with a seeded random generator.
3. Count the runs where the event happens.
4. Estimate p̂ = count / N.
5. Report the error: SE = sqrt(p̂ (1 - p̂) / N).

Step 5 is what separates a simulation from a guess. By the CLT (previous lesson), p̂ is roughly normal around the true p, so p̂ ± 2 SE is an approximate 95% range.

Worked example: the birthday problem

In a room of 23 people, what is the probability that at least two share a birthday (ignoring leap years and assuming all 365 days are equally likely)? The exact answer is about 0.507, which surprises most people. Here is a simulation using only the standard library:

import random

def shared_birthday(n_people, rng):
    days = [rng.randrange(365) for _ in range(n_people)]
    return len(set(days)) < n_people       # a repeat shrinks the set

rng = random.Random(2026)                  # seeded: reproducible
trials = 20_000
hits = sum(shared_birthday(23, rng) for _ in range(trials))
p_hat = hits / trials
se = (p_hat * (1 - p_hat) / trials) ** 0.5
print(p_hat, se)                           # about 0.507, SE about 0.0035

With 20,000 trials the standard error is about 0.0035, so the estimate is reliable to roughly ±0.007. You can check it against the exact formula:

P(no shared birthday) = (365/365) × (364/365) × ... × (343/365) ≈ 0.493
P(shared)             ≈ 1 - 0.493 = 0.507

Worked example: estimating π

Throw random points into the unit square. The share that lands inside the quarter circle x² + y² ≤ 1 is π/4:

rng = random.Random(1)
N = 100_000
inside = sum(rng.random() ** 2 + rng.random() ** 2 <= 1 for _ in range(N))
print(4 * inside / N)       # about 3.14

This is a toy, but the same idea (estimate an area or integral by random sampling) is how physicists, engineers and Bayesian statisticians compute quantities that have no closed form.

Worked example: tolerance stacking

An assembly is three parts end to end, each Normal(mean 10.0 mm, SD 0.05 mm), and it must be at most 30.12 mm long. With NumPy:

import numpy as np
rng = np.random.default_rng(0)
total = rng.normal(10.0, 0.05, size=(100_000, 3)).sum(axis=1)
print((total > 30.12).mean())      # about 0.083

Here a formula exists (the sum is normal with SD 0.05 × sqrt(3) ≈ 0.087), and the simulation agrees with it. But if one part were uniform, or the parts were correlated, the simulation would need only a one-line change while the formula would need new math.

Choosing N

The error falls as 1/sqrt(N). For a probability near 0.5:

Rare events need more runs. To pin down a 1-in-1,000 probability to within about ±10% of its value (a 95% range), you need around 400,000 runs. Analysts use variance-reduction tricks (importance sampling, stratified sampling) for very rare events.

Common mistakes

• No seed. Results change on every run and nobody can reproduce your figure.
• No error bar. "The probability is 0.507" from 200 trials (SE ≈ 0.035) is overconfident.
• Modelling the wrong process. A simulation is only as good as its assumptions. Birthdays are not quite uniform across the year; say so.
• Reusing the same random numbers by accident, such as reseeding inside the loop, which makes every trial identical.

Recap

• Monte Carlo: model one run, repeat N times, count, divide.
• Always seed the generator and report SE = sqrt(p̂(1 - p̂)/N).
• Error shrinks with 1/sqrt(N); rare events need many runs.
• Use simulation to answer questions with no formula and to check the ones that have a formula.

Try it

In the Do tab, simulate the birthday problem and compare your estimate with the exact answer.