Monte Carlo Simulation

In the 1940s, working on the atomic bomb, Stanislaw Ulam fell ill and passed the time playing Canfield solitaire. He wondered what fraction of games could actually be won — and found the combinatorics hopeless. Then it struck him: rather than calculate the odds, why not simply deal a hundred hands and count how many came out? His colleague John von Neumann seized on the idea for the fearsome neutron-diffusion integrals nobody could do by hand, and, needing a secret code name, someone suggested Monte Carlo after the casino where Ulam's uncle gambled. The joke stuck, and one of the most widely used numerical methods in science was born.

The idea is almost embarrassingly simple. To estimate a quantity you cannot compute directly, draw many random samples and average the results. The average is your estimate; the more samples you draw, the closer it creeps to the truth. That is the whole method — and it works on problems from pricing options to sizing hospitals to rendering film.

Why averaging random samples works

Suppose the answer you want is the expected value of some random quantity X — a mean, a probability, an integral, a risk. You can't evaluate \mathbb{E}[X] in closed form, but you can generate independent samples X_1, X_2, \dots, X_n and take their sample mean:

\hat{\theta}_n = \frac{1}{n}\sum_{i=1}^{n} X_i \;\xrightarrow[n\to\infty]{}\; \mathbb{E}[X].

The law of large numbers guarantees that this sample average converges to the true expectation. The central limit theorem goes further and tells you how fast: the standard error of the estimate shrinks like

\text{error} \;\approx\; \frac{\sigma}{\sqrt{n}},

where \sigma is the standard deviation of a single sample. That 1/\sqrt{n} is the heartbeat of Monte Carlo — and, as we'll see, both its blessing and its curse.

The canonical example: estimating π by throwing darts

Draw a square of side 2 centred on the origin, and inscribe a circle of radius 1. The square has area 4; the circle has area \pi. Now throw darts uniformly at random into the square. The chance a dart lands inside the circle is the ratio of areas, \pi/4. So if a fraction p of your darts land inside, then \pi \approx 4p.

No calculus, no integral — just counting. Throw 1000 darts, find that 786 land inside, and you estimate \pi \approx 4 \times 0.786 = 3.144. Every quantity you cannot compute directly, from the volume of a strange region to the value of a financial derivative, can be recast as "what fraction of random trials succeed?" or "what is the average payoff of a random draw?"

Watch it converge

Here is a Monte Carlo estimate of \pi as the number of samples n grows. Each dot is the running estimate after that many darts; the shaded band is the \pm c/\sqrt{n} error envelope. The estimate rattles around early on and then funnels toward the true value \pi = 3.14159\ldots:

Notice the shape of the funnel. To go from a handful of correct digits to one more digit, you don't need a bit more work — you need the whole envelope to shrink again, and because it shrinks like 1/\sqrt{n}, that costs a fourfold increase in samples every time.

The same trick estimates integrals and risks

Any integral is an average in disguise. To estimate I = \int_a^b g(x)\,dx, draw points x_1,\dots,x_n uniformly on [a,b] and average:

I \;\approx\; (b-a)\,\frac{1}{n}\sum_{i=1}^{n} g(x_i).

In high dimensions this is transformative. A grid of 100 points per axis needs 100^{d} evaluations in d dimensions — hopeless past a handful of variables — yet Monte Carlo's 1/\sqrt{n} error does not depend on the dimension at all. That is why it dominates in finance (a portfolio with hundreds of risk factors) and physics (the many-particle integrals that computational physics was built to attack). To estimate a risk — say, the chance a project overruns its budget — you simulate thousands of possible futures with random costs and durations, and count the fraction that blow the budget.

Yes — cleverly. The raw error is \sigma/\sqrt{n}, and you can attack either factor. Throwing more darts shrinks \sqrt{n}, but slowly. The smarter route is variance reduction: shrink the \sigma instead. Antithetic variates pair each random draw u with its mirror 1-u; the two are negatively correlated, so their average has less variance than two independent draws. Control variates lean on a similar quantity whose answer you already know, importance sampling concentrates draws where the integrand is big, and stratified sampling forces even coverage. A good variance-reduction scheme can do the work of ten or a hundred times as many raw samples — the same accuracy for a fraction of the compute.

Because the error scales like 1/\sqrt{n}, to halve it you must quadruple the number of samples; to get one extra decimal digit (ten times smaller error) you need 100\times the samples. Monte Carlo is wonderfully robust — it barely cares about dimension, smoothness or the shape of the region — but it is slow to high precision. If you need six digits, a deterministic method (quadrature, a series) will usually crush it. Reach for Monte Carlo when the problem is high-dimensional, awkwardly shaped, or genuinely random — not when you simply want many digits cheaply.