Some probability problems have an elegant solution on paper. Many others do not: too many variables, complicated dependencies, irregular distributions. For these there is a method of almost disarming simplicity: simulate the phenomenon thousands of times with random numbers and see what happens. It is the Monte Carlo method, and with a language like Python a few lines are enough to use it.
The idea
If we cannot compute a probability or an expected value, we can:
- write a program that generates one possible scenario according to the rules of the problem;
- repeat it a great many times with different random numbers;
- count how often the event we care about happens, or compute the average of the result.
By the law of large numbers, the observed frequency approaches the true probability and the observed average the expected value. The method works for any problem that can be described as a random process, even when no formula is available.
First example: estimating π
A classic example, useful for understanding the mechanism, is estimating π. Uniform random points are generated in the unit square. The fraction falling inside the quarter circle of radius 1 tends to the ratio of the areas, that is π/4.
import numpy as np
rng = np.random.default_rng(42)
n = 1_000_000
x, y = rng.random(n), rng.random(n)
estimate = 4 * np.mean(x**2 + y**2 <= 1)
print(estimate)
With a million points the estimate is usually correct to the second or third decimal place. It is certainly not the best way to compute π, but it shows well how a proportion of random events can estimate a deterministic quantity.
Second example: the duration of a project
Here is a case more useful in real work. A project consists of three tasks in sequence, and for each we have an estimate with a minimum, a most likely value and a maximum, in days:
| Task | Minimum | Most likely | Maximum |
|---|---|---|---|
| Analysis | 2 | 3 | 5 |
| Development | 4 | 6 | 10 |
| Testing | 1 | 2 | 4 |
Adding up the most likely values gives 11 days, and it is easy to promise that date. Instead, let us simulate the project by modelling each task with a triangular distribution:
n = 100_000
duration = (rng.triangular(2, 3, 5, n)
+ rng.triangular(4, 6, 10, n)
+ rng.triangular(1, 2, 4, n))
print(duration.mean()) # about 12.3 days
print(np.mean(duration <= 11)) # about 0.2
print(np.percentile(duration, [50, 90])) # about 12.3 and 14.4
The result is instructive: the probability of finishing within 11 days is only about 20%. The reason is that the estimates are asymmetric: each task can run late by much more than it can finish early, and delays add up. To have a 90% probability of meeting the deadline you would have to quote about 14.4 days. A few lines of simulation turn an optimistic point estimate into a realistic distribution of outcomes.
How precise a simulation is
The result of a Monte Carlo simulation is itself an estimate, with its own error. For a probability p estimated with n trials, the standard error is about √(p(1 − p)/n). For a mean, it is the standard deviation of the result divided by √n.
The consequence is the same as for sampling: to halve the error you need four times as many simulations. With 10,000 trials a probability around 20% is estimated with a standard error of about 0.4 percentage points; with a million trials, about 0.04.
It is always worth reporting this error, or at least repeating the simulation with different seeds to see how much the results fluctuate.
Good practice
Use a modern generator and a seed
In NumPy the recommended way is to create a generator with np.random.default_rng(seed). Fixing the seed makes the simulation reproducible: anyone rerunning the code gets the same numbers, which is essential for checking and correcting an analysis.
Vectorise
Generating all random numbers at once with NumPy arrays is much faster than using Python loops. In the examples above, a hundred thousand scenarios are computed in a fraction of a second.
Validate on a known case
Before trusting a complex simulation, test it on a simplified version whose exact solution you know. If the simulation of rolling two dice does not return 7 as the average sum, there is a bug in the code.
Typical mistakes
- Ignoring dependencies: if two tasks tend to run late together, because they depend on the same person or the same supplier, simulating them as independent underestimates the risk.
- Choosing distributions arbitrarily: the result cannot be better than the input assumptions. If the minimum and maximum estimates are too optimistic, the simulation will be optimistic with great precision.
- Confusing precision with accuracy: a million simulations reduce random error, not the error due to a wrong model.
- Using unsuitable generators where security matters: pseudo-random generators for simulation are not suitable for creating keys or tokens; for that, cryptographic generators are needed.
In short
- Monte Carlo simulation estimates probabilities and expected values by repeating a random process many times.
- It works even when no formula is available, as long as the process can be described.
- Adding up the most likely values of asymmetric estimates leads to optimistic forecasts.
- Simulation error decreases with √n: halving it requires four times as many trials.
- Fixed seeds, vectorised code and validation on known cases make the simulation reliable.
