Monte Carlo methods: answering hard questions by rolling dice

There is a class of question that has a clean mathematical answer you will never write down. What is the chance this project finishes before March? What range should the revenue forecast carry, given that three of its inputs are themselves guesses? What is the area of this awkward shape?

Monte Carlo is the move you make when the pen fails. Instead of solving for the answer, simulate the situation many times, watch what happens, and count. The average of the simulations is your estimate. The spread of the simulations is your uncertainty, free of charge.

It sounds like giving up. It is closer to the opposite: it trades a problem you cannot solve for one a computer can grind through, and it comes with an honest, predictable error. The method is old, it is boring, and it quietly underpins a large share of applied statistics.

The toy case: estimating pi

Draw a unit square. Draw a quarter circle inside it. The square has area 1, the quarter circle has area pi/4. Now throw darts at the square at random. The fraction that land inside the arc estimates pi/4, so four times that fraction estimates pi.

import numpy as np

rng = np.random.default_rng(0)
n = 1_000_000
x, y = rng.random(n), rng.random(n)
inside = (x**2 + y**2) <= 1.0

pi_hat = 4 * inside.mean()
stderr = 4 * inside.std(ddof=1) / np.sqrt(n)
print(f"{pi_hat:.4f} +/- {1.96 * stderr:.4f}")

A million darts gets you about 3.14 give or take 0.002. Nobody computes pi this way, and that is the point of the example: because the true answer is known, you can see the method working and see exactly how fast.

Both halves of that are worth seeing:

Left: 300 darts, 234 of them inside the arc, which gives 3.12 rather than 3.1416. Right: what it costs to do better. The typical error is 1.64 over the square root of n, so going from 625 darts to 2,500 halves it, and 2,500 to 10,000 halves it again. Nothing about the problem changes that rate. The structure generalises. Every Monte Carlo estimate is an average of a function evaluated at random points. Whether the function returns “did the dart land inside” or “did the project ship on time” or “what was the profit in this scenario”, the machinery is identical: sample, evaluate, average.

Why the error falls as one over root n

Each dart is an independent draw. The average of n independent draws has a standard error of sigma over the square root of n. That single fact is the whole cost model of Monte Carlo, and it has two consequences worth memorising.

Accuracy is expensive. To halve the error you need four times the samples. For one more decimal place — ten times more accurate — you need a hundred times more samples. Anyone promising six-digit precision from a simulation is either running an enormous job or not measuring their error.

Dimension is free. The one-over-root-n rate does not care whether you are integrating over 2 variables or 200. Grid methods do care, badly: a grid with 10 points per axis needs 10^200 cells in 200 dimensions. This is the reason Monte Carlo dominates in high dimensions rather than being a last resort. In one dimension, use calculus or a quadrature rule. In twenty, sample.

The practical reading of both points together: pick the accuracy you actually need before choosing n. A go/no-go decision that turns on whether risk is above or below 5% does not need the third decimal place. Ten thousand runs usually settles it; ten million is often someone soothing themselves.

What it is really for

The dart-throwing framing undersells it. Three uses carry most of the weight in practice.

Propagating uncertainty. You have a model with uncertain inputs. Rather than guessing what the output uncertainty must be, push distributions through the model and look at the output distribution. Five lines:

price = rng.normal(50, 4, 100_000)
units = rng.normal(2_000, 300, 100_000)
cost = rng.normal(60_000, 5_000, 100_000)
profit = price * units - cost
print(np.percentile(profit, [5, 50, 95]), (profit < 0).mean())

That last number — the probability of a loss — is the thing a decision maker wants and cannot get from a spreadsheet of point estimates. This is also far better than the usual “best case / worst case / base case” three-column theatre, which quietly assumes every input goes wrong at once.

Answering questions about a statistic you cannot derive. The bootstrap is Monte Carlo: resample the data, recompute the statistic, look at the spread. Permutation tests are Monte Carlo. So is power analysis for an experiment design that has no textbook formula — simulate the effect you expect, run your own analysis on the simulated data, count how often you detect it.

Simulating a process over time. Queues, inventory, cash runway, retries under load. Anything where the state at each step depends on the last and the rules are easy to write but hard to solve.

The failure mode

Monte Carlo answers the model you gave it, with impressive precision, whether or not the model is any good. Run enough samples and you will produce a tight interval around a number that is wrong.

Two ways this bites. First, the input distributions are usually guesses wearing a normal curve. If revenue has never once been symmetric around its mean, sampling it from a normal understates the bad tail — which is the tail anyone cares about. Second, and more common: sampling each input independently when the inputs move together. In the profit example above, price and units are almost certainly correlated in the real world, and pretending otherwise makes the simulated range look narrower and friendlier than reality.

So sanity-check the output the way you would any model. Does the simulated distribution look like history? Do the extreme runs describe situations that could actually happen? Fix a seed so results are reproducible, and re-run with a different one to confirm the answer is stable. Report the interval, never the mean alone — and when you present it, say plainly what the simulation assumed, because that is the part stakeholders can challenge.

Monte Carlo Simulation by MarbleScience is a short visual introduction if the sampling idea still feels abstract.

One caveat about the whole method. It is only as good as the distributions fed into it, and for a genuinely new product there is no history to fit them to. Sampling from invented distributions produces an interval that looks like evidence and is entirely assumption, which is a worse position than a single guessed number, because at least a guessed number is read as a guess. Where there is no data, say so, put the range on the assumption instead, and argue about that.