Monte Carlo Methods: Getting the Right Answer by Guessing

Answering a question by throwing dice at it

A Monte Carlo method replaces a calculation you cannot do with an experiment you can repeat. Write the quantity you want as an expected value, draw random samples from the right distribution, and take the average. The law of large numbers guarantees the average converges to the true value; the central limit theorem tells you how fast.

The idea comes from Los Alamos in the 1940s. Stanislaw Ulam, recovering from illness, tried to work out the probability that a game of solitaire comes out — and realised it was far easier to deal a hundred hands and count than to enumerate the combinatorics. John von Neumann saw immediately that the same trick applied to neutron transport, and Nicholas Metropolis supplied the code name, after the casino in Monaco where Ulam's uncle liked to gamble.

π from a square and a circle

The canonical demonstration, and the one running in the simulation on this site, is estimating π. Take the unit square and the quarter circle inscribed in it. Drop points uniformly at random into the square. The probability that a point lands inside the quarter circle is the ratio of the areas, which is π/4. So count the hits, divide by the throws, multiply by four.

let inside = 0;
for (let i = 0; i < n; i++) {
  const x = Math.random(), y = Math.random();
  if (x * x + y * y <= 1) inside++;      // no sqrt needed
}
const piEstimate = 4 * inside / n;

Nothing about this is specific to circles. Any integral is an area, and any area is the probability that a uniform point lands under the curve — so this same three-line loop is a general-purpose integrator. The only thing that changes from problem to problem is what you sample and what you count.

The 1/√N law, and why it is both terrible and wonderful

Each throw is a Bernoulli trial with success probability p = π/4 ≈ 0.7854, so its variance is p(1 − p) ≈ 0.1685. The estimator is an average of N such trials, so its standard error is √(p(1−p)/N), and the estimate of π — four times that — has a typical error of

σ(π̂) = 4 · sqrt(p(1 − p) / N) ≈ 1.64 / sqrt(N)

N = 10⁴    → typical error ≈ 0.016     (about 2 correct digits)
N = 10⁶    → typical error ≈ 0.0016    (about 3 correct digits)
N = 10⁸    → typical error ≈ 0.00016   (about 4 correct digits)

That is the defining property of Monte Carlo: the error falls as 1/√N, so buying one more decimal digit costs a hundred times more samples. As a way of computing π it is hopeless — Machin-like series get more digits in a handful of terms than this loop gets in a billion throws.

But look at what is not in that formula: the dimension. A deterministic quadrature rule on a grid needs mᵈ points to put m points along each of d axes, and its accuracy collapses as d grows — the curse of dimensionality. Monte Carlo's 1/√N does not depend on d at all. In two dimensions it is a bad joke; in the twenty or two hundred dimensions of a light-transport integral, a financial payoff or a statistical-mechanics partition function, it is often the only thing that works. Every path tracer, every option pricer and every Bayesian posterior sampler is built on that one observation.

Variance reduction: the only lever you have

Since the exponent in 1/√N is fixed, the only way to get a better answer for the same cost is to shrink the constant in front — the variance of what you are averaging. Four standard techniques:

importance sampling   sample where the integrand is large, then divide
                      by the sampling density (the weight). The classic
                      case: sampling a light source directly instead of
                      firing rays into a mostly black hemisphere.

stratified sampling   split the domain into K strata and sample each,
                      forbidding the clumps and gaps of pure randomness.

antithetic variates   for every sample u, also use 1 − u. The negative
                      correlation cancels part of the error.

control variates      subtract a correlated quantity whose exact mean
                      you already know, and add that mean back.

There is also quasi-Monte Carlo, which abandons randomness entirely for a low-discrepancy sequence (Sobol, Halton) — deterministic points engineered to fill space more evenly than random ones ever do. For smooth, moderate-dimensional integrands its error can approach 1/N rather than 1/√N, which is an enormous win; the price is that the classical error bars no longer apply, and the advantage erodes as the dimension climbs. Randomised QMC (scrambling the sequence) buys the error estimate back.

When you cannot sample the distribution directly

Everything above assumes you can draw independent samples from the target distribution. In Bayesian statistics and statistical physics you usually cannot: you know the density only up to an unknown normalising constant. The answer is Markov chain Monte Carlo — build a random walk whose long-run stationary distribution is the target, then average along the walk. The Metropolis–Hastings rule is startlingly simple: propose a move, and accept it with probability

α = min(1, π(x_new) / π(x_old))       // symmetric proposal

accept  → the chain moves
reject  → the chain stays where it is, and you count x_old AGAIN

The normalising constant cancels in the ratio, which is the entire point. The cost is that consecutive samples are correlated, so N MCMC steps carry less information than N independent draws, and you must discard an initial burn-in while the chain forgets where it started. Gibbs sampling, Hamiltonian Monte Carlo and the No-U-Turn Sampler are refinements of the same skeleton.

The random numbers themselves

A Monte Carlo result is only as good as its generator. Three practical rules. First, use a well-tested PRNG — a modern counter-based or xorshift-family generator, not a home-made linear congruential one, whose low bits are notoriously non-random and whose points famously fall in hyperplanes. Second, seed it explicitly and record the seed: a Monte Carlo run that cannot be reproduced cannot be debugged. Third, watch the period and the correlation between streams if you parallelise — two workers accidentally sharing a stream is a silent, catastrophic bug that no test will catch, because the answer still looks plausible.

Math.random() in a browser is fine for a visualisation like this one (it is not seedable and not cryptographically strong, but it is fast and its statistical quality is adequate); for a scientific run, use a generator you can seed and reproduce.

Frequently asked questions

Try it live

Everything above runs in your browser — open Monte Carlo π and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.