Learn AI Create the Future — with Pranay
Main site

Resampling, Likelihood and Bayes

Monte Carlo simulation

When the formula is too hard, roll the dice instead — simulate the process thousands of times and count what happens.

Read these first

On this page 5
  1. Why it exists
  2. How it works
  3. A real example you have seen
  4. Remember this
  5. What to learn next

One lesson, three depths. Pick the one that fits you today — you can switch any time.

Beginner — No maths. Plain English.

Monte Carlo simulation answers hard probability questions by acting the situation out thousands of times and counting the results.

How do you find out if a coin is bent? Flip it. A thousand flips beats any argument about its shape.

Now the leap: any question about chance can be "flipped". Cannot work out the odds on paper? Build the situation inside a computer, run it ten thousand times, and count. The counting is the answer.

Why it exists

During the 1940s, physicists at Los Alamos needed to know how neutrons travel through material. The exact equations were hopeless. Stanislaw Ulam, recovering from illness, played solitaire and wondered about his odds of winning. He could not compute them. But he could deal a hundred games and count.

He and John von Neumann turned that into a method. The code name honoured the Monte Carlo casino, where Ulam's uncle gambled. The method now runs inside banks, weather models, and game engines.

Its superpower is handling tangled randomness. One uncertain thing has formulas. Five uncertain things feeding into each other — that is where paper fails and simulation shines.

How it works

describe each uncertainty   →  "task A takes 2 to 6 days, likely 3"
wire them together          →  "B starts after A; C runs alongside"
run the whole story once    →  one possible future:  6.8 days
run it 100,000 times        →  100,000 possible futures
count                       →  "finishes within 7 days in 61% of them"

Each run draws one random value per uncertainty and plays the story out. No single run means anything. The pile of runs is a map of every way things could go, weighted by how likely each way is.

The answer's precision grows with more runs — slowly. Ten thousand runs give roughly one extra digit over one hundred. Patience is the price of generality.

A real example you have seen

Weather apps saying "70% chance of rain" run dozens of atmosphere simulations with slightly different starting guesses, and count how many turn rainy. Cricket win-probability tickers simulate the remaining overs thousands of times. Pension planners simulate market futures to ask "will the savings last?".

Remember this

  • Cannot solve it? Simulate it and count.
  • One run is a story; ten thousand runs are a map of possible futures.
  • Precision grows slowly with more runs — generality is the real payoff.

What to learn next

Developer — Code and libraries.

Setup

bash
pip install numpy

Outputs verified with numpy 1.26, CPU. Digits reproduce with this version; other versions may shift the last decimal.

Will the project finish on time?

Three tasks with uncertain durations. Tasks 1 and 2 run one after the other; task 3 runs alongside them. The deadline is 7 days. No clean formula exists for this — which is the point.

project_mc.py
import numpy as np

rng = np.random.default_rng(11)
n = 100_000
t1 = rng.triangular(2, 3, 6, n)      # best 2, most likely 3, worst 6
t2 = rng.triangular(1, 2, 4, n)
t3 = rng.triangular(3, 5, 10, n)
total = np.maximum(t1 + t2, t3)      # done when the slower path is done

print(f"mean finish time: {total.mean():.1f} days")
print(f"P(done within 7 days) = {(total <= 7).mean():.2f}")
print(f"90% of runs finish within {np.quantile(total, 0.9):.1f} days")
Output
mean finish time: 6.7 days
P(done within 7 days) = 0.61
90% of runs finish within 8.3 days

Fifteen lines, and you have answers a spreadsheet full of "best case / worst case" columns cannot give.

The walkthrough

rng.triangular(2, 3, 6, n) draws 100,000 possible durations for task 1 at once — no Python loop. The triangular distribution is the workhorse of quick modelling: ask a colleague for best, most-likely and worst, and you have its three parameters.

np.maximum(t1 + t2, t3) is where the tangling lives. The project ends when the slower of the two paths ends. This single line is what breaks pen-and-paper: the maximum of random quantities has no friendly algebra, but the simulation does not care.

Notice the trap the mean sets. Adding the most-likely values gives 3 + 2 = 5 days on one path, 5 on the other — "about 5 days, well inside 7". The simulation says the deadline fails 39% of the time. Randomness piles up in ways point estimates hide; this gap is called the flaw of averages.

(total <= 7).mean() is the counting step: a boolean array's mean is the fraction of True — the estimated probability in one idiom.

How precise is 0.61? The standard error of a counted probability is about sqrt(p*(1-p)/n) — here 0.0015. Report 0.61 and move on; chasing a third digit needs a hundred times more runs.

Common mistakes

Reporting only the mean of the pile. The decision-grade numbers are the probabilities and quantiles — "61% on time", "90% done by 8.3 days". The mean alone re-creates the flaw of averages you came here to escape.

Simulating dependent things as independent. If rain delays task 1, it likely delays task 3 too. Independent draws understate risk — the same mistake that broke mortgage models in 2008. Model shared causes explicitly: draw a "weather" variable first, and let both tasks depend on it.

Reusing one seed everywhere and calling the result "the answer". Run with two or three seeds once; if answers move more than you can tolerate, raise n. Then fix one seed for reproducibility — see random seeds and reproducibility.

Looping in Python when arrays would do. A for loop over 100,000 runs is a hundred times slower than the vectorised draw above. Shape: one array per uncertainty, length n.

Try it yourself

Add a fourth task: after everything else, deployment takes between 0.5 and 2 days, uniformly (rng.uniform). Rerun and find the new on-time probability — then find what deadline you could promise with 90% confidence.

What to learn next

Researcher — Mathematics and papers.

The estimator and its guarantees

Monte Carlo estimates $\mu = \mathbb{E}[f(X)]$ by

$$ \hat\mu_n = \frac{1}{n} \sum_{i=1}^{n} f(X_i), \qquad X_i \sim p \ \text{i.i.d.} $$

Where:

  • $X$ — the random inputs with density $p$; $f$ — the quantity computed per run.
  • Probabilities are the special case $f = \mathbf{1}[\text{event}]$.

The law of large numbers gives consistency; the CLT gives the error law:

$$ \hat\mu_n - \mu \;\approx\; \mathcal{N}!\left(0, \frac{\sigma^2}{n}\right), \qquad \sigma^2 = \operatorname{Var}(f(X)) $$

Convergence is $O(n^{-1/2})$ independent of dimension — the property that makes Monte Carlo the only game in town for high-dimensional integrals, where deterministic quadrature's error grows exponentially with dimension. The price: each extra digit of precision costs 100× the samples.

Variance reduction

The constant $\sigma$ is attackable even though the rate is not:

  • Antithetic variates: pair each draw $U$ with $1-U$; negative correlation between pairs cancels variance for monotone $f$.
  • Control variates: with a correlated $g$ whose mean is known, estimate $\mathbb{E}[f] - \beta(\mathbb{E}[g] - \bar g)$; optimal $\beta = \operatorname{Cov}(f,g)/\operatorname{Var}(g)$.
  • Importance sampling: draw from a proposal $q$ concentrated where $f p$ is large and reweight by $w = p/q$:

$$ \mathbb{E}_p[f(X)] = \mathbb{E}_q!\left[f(X) \frac{p(X)}{q(X)}\right] $$

essential for rare events, where naive counting sees almost no hits; a poor $q$ can make variance infinite, so diagnose with the effective sample size $\left(\sum w_i\right)^2 / \sum w_i^2$.

  • Stratified sampling and its scaled-up form, quasi-Monte Carlo: low-discrepancy sequences (Sobol', Halton) achieve near-$O(n^{-1})$ error for smooth integrands in moderate dimension (scipy.stats.qmc).

Random number foundations

Modern simulation rests on counter-based and permuted-congruential generators; numpy's default PCG64 has period $2^{128}$ and passes TestU01's BigCrush. Reproducibility contract: np.random.default_rng(seed) produces identical streams across platforms for a fixed numpy version. Parallel workers need rng.spawn() or SeedSequence — never worker-index seeds like 1, 2, 3, which can correlate streams.

Where this sits in ML

Monte Carlo is the substrate of much of the field: stochastic gradient descent estimates the full-data gradient by minibatch sampling; dropout at inference (MC dropout, Gal and Ghahramani, 2016) estimates predictive uncertainty; policy-gradient reinforcement learning estimates expected return by rollout; MCMC extends the idea to distributions you can only score, not sample.

History

Metropolis and Ulam (1949), The Monte Carlo method, Journal of the American Statistical Association — the public debut. Buffon's needle (1777) is the recognised ancestor: estimating $\pi$ by dropping needles and counting line crossings.

What to learn next