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.