The Monte Carlo method

How many square kilometres of water are there in the Kapchagay reservoir? There is no need to walk the shore: take a satellite image and throw points onto it at random. The share of points that landed in water, times the area of the frame, is the answer: approximate, but with an honestly stated error. The frame in the scene is 3517.5 km², and with each new point the band around the estimate on the chart narrows.

Look at the chart. After a thousand points the estimate is 1097 km² plus or minus a hundred. After a million, 1124 ± 3. A thousand times more points, and the error shrank thirtyfold: √1000 ≈ 31.6. This is the main law of the method, and the whole article is about it. The second thing on the chart is the dashed “handbook” line at 1847 km². A million points do not approach it and never will. Why is in the section on the limits of the method; press “What counts as water” and look at the eastern end of the reservoir.

The idea

If the number sought is a probability or the mean of something random, it can be estimated by experiment: play out the chance n times and take the share of successes or the mean of the results. The error of such an estimate falls as σ/√n, where σ is the spread of one trial. The formula contains neither the dimension of the problem nor the complexity of the shape: only the spread and the number of trials.

The area of the reservoir is the mean of the function “1 if the point is in water, 0 otherwise” over the whole frame, times the area of the frame. In exactly the same way the integral of any function over a region is the mean of its values at random points times the volume of the region. The method needs neither a formula for the boundary nor its smoothness: only the ability to answer whether a point is inside or not. For an image that is the colour of a pixel, for a complex part a check against the drawing, for a neutron in a reactor a playing out of its fate.

Why 1/√n

For a point “in water or not”, σ² = p(1 − p), where p is the true share of water, about 0.32 here. By the central limit theorem, for large n the estimate is almost normally distributed, and in 95% of experiments it deviates from the truth by at most 1.96 · σ/√n. Hence the familiar interval p̂ ± 1.96 · √(p̂(1 − p̂)/n). To narrow it by half takes four times as many points; to get one more correct decimal digit, a hundred times as many.

Slow? Yes. But the law does not depend on how many dimensions the problem lives in. In two dimensions the same accuracy can be had with a grid, √n by √n points. But in thirty, a grid with two points along each axis is already a billion points, and one point along an axis is no grid at all. Random points do not care: σ/√n.

Is the interval honest

“95%” promises this: if the whole experiment is repeated again and again, the interval will cover the truth in 95 cases out of a hundred. Let us check the promise directly: a hundred independent experiments of n points each, each with its own interval. The red ones missed the truth.

For the reservoir at n = 100 the familiar formula keeps its word: about five misses per hundred. Now choose “pond, 1%”, a body of water that takes up a hundredth of the frame, and n = 30. Three experiments out of four will not catch a single point in water, and the formula will give an interval of zero width: “there is no water, and we are sure of it”. Exactly how often can be computed without any simulation: the number of hits follows the binomial law, and the probability of covering the truth is the sum of the probabilities of those k whose intervals cover it.

   share      n    P(no hit)   textbook   Wilson
    lake    100        0.000      0.946    0.959
    pond     30        0.740      0.260    0.964
    pond    100        0.366      0.633    0.921

95% is promised, but for the pond with 30 points the familiar interval covers the truth in 26% of cases. In 2001 Lawrence Brown, Tony Cai and Anirban DasGupta showed that the coverage of this formula jumps erratically even for large n, and that textbook advice on when it can be trusted is unreliable, and they recommended the Wilson interval instead.

Edwin Wilson’s idea (1927) is not to substitute the estimate p̂ into the formula for the spread but to ask honestly: for which values of p is the observed share p̂ plausible, that is, differs from p by at most 1.96 spreads computed from p itself?

Square it, and you get a quadratic inequality in p, whose solution is an interval with ends

With p̂ = 0 it does not collapse to a point, never goes outside [0, 1] and keeps coverage close to what is promised: for the same pond, 96.4%. Switch the scene to “Wilson” and compare.

Code

The program reads the water mask, the same image reprojected to an equal-area projection where each pixel covers the same 97 × 97 m, and throws points onto it. The random number generator is our own and the same in Python and Go, so both versions throw the same points; the first rain in the scene is them too.

frame 960x390 px = 3517.5 km2; water pixels 119639 = 1124.0 km2
        n   estimate, km2   ±95%      Wilson 95%         covers
       10          1055.2    999.1   [ 379.2, 2121.8]   yes
     1000          1097.4    101.0   [ 999.1, 1200.8]   yes
  1000000          1124.2      3.2   [1120.9, 1127.4]   yes

Ten points give “1055 ± 999”: honest but useless. The Wilson interval on the same ten points, from 379 to 2122, is asymmetric: it remembers that an area cannot be negative.

What it costs

The time is O(n) trials, each costing one “inside or not” check, and the memory O(1): only a counter needs keeping. Accuracy is what is paid for: the error is ∝ 1/√n, and each next digit costs a hundred times more. A million points take seconds even in Python; two more digits would take ten billion.

This law can be got round by giving up full independence of the points. Cut the frame into k × k cells and throw one random point into each. The error now comes only from the cells the shore passes through: cells entirely in water or entirely on land give the exact answer. There are about k such “shore” cells, that is, √n, and the error falls as n^−3/4 instead of n^−1/2; this result for pixels with an edge is in Don Mitchell’s work (1996), and the technique itself was brought into computer graphics by Robert Cook, Thomas Porter and Loren Carpenter in 1984. Switch the first scene to “By cells”: with the same n the interval is noticeably narrower.

     n   plain: rms error, km2   stratified   ratio
   100                   139.7         85.4     1.6
  6400                    22.7          4.1     5.6

On 6400 points the cells err 5.6 times less, and the gap grows with n. In many-dimensional problems cells are hard to build, and their relatives are used there: quasi-random sequences, which lay points out evenly without a grid.

Where it cannot be done without

An area on a plane can be computed without randomness too, with a grid or by the polygon formula. The Monte Carlo method is irreplaceable where there are many dimensions. Here is a problem with ninety dimensions. A million tenge is kept in dollars, euros and yuan in equal parts. How much will it be worth in tenge in a month? A month is 30 days, and each day three rates change: the outcome depends on ninety random numbers. There is no formula: daily changes in rates obey no tidy law, and a grid with two points on each of 90 axes is 2⁹⁰ scenarios.

Monte Carlo acts simply. There are five years of official rates, 1899 daily changes. Assemble a month from 30 randomly chosen days of history, each day whole, all three currencies at once, so that the dollar and euro move together as they really did, and compute what the piggy bank is now worth. Repeat twenty thousand times. Such a technique, where the experiment is played out on reshuffled real data, is called the bootstrap.

1899 daily moves, 20000 simulated months
mean after a month: 1,000,933 tenge
P(lose more than 5%) = 3.95% ± 0.27%
5% worst month: below 953,139 tenge, 95% interval [952,352, 954,076]
real 30-day windows losing more than 5%: 87 of 1870 = 4.65%

The probability of losing more than 5% in a month is 3.95% ± 0.27%. The worst 5% of months leave less than 953 thousand of the million; this bound has an interval too: the outcomes neighbouring in rank, from 952 to 954 thousand. This is how banks estimate the market risk of a portfolio: they play out scenarios on history.

And now the last line: in real history months with such a loss make up 4.65%, more than the simulation counted, and outside its interval. Monte Carlo made no arithmetic mistake: it honestly answered a question about a model in which days are independent. The main reason, apparently, is that real bad months are runs of bad days in a row, as in the spring of 2022, when within a month the dollar rose by 20% and then fell by 12%, while reshuffling days breaks such runs up. The interval describes random error and knows nothing about whether the model is right.

Where the method is powerless

Systematic error. Back to the reservoir. A million points gave 1124 ± 3 km², while the handbook says 1847 km². The method is not to blame: it converges to what the classifier counts as water, 1124.0 km² over all pixels of the image. A dark blue-green pixel counts as water, and the murky shallow eastern end, where the river flows in, looks like land on the 2016 image. And which water level the handbook’s 1847 km² refers to, the handbook does not say; it is known that the reservoir was never filled to its design volume. No number of added points will reduce such an error: it lies in the definition, not in the randomness.

Rare events. To estimate a probability of 10⁻⁶ with a relative error of 10%, the event must happen at least a hundred times: a hundred million trials. The relative error of an estimate of a small probability p is about 1/√(np): for p = 10⁻¹² the method is powerless. Such probabilities are computed exactly, as in the second problem, or the distribution is changed so that the rare becomes frequent, and the weights are recomputed: this is importance sampling.

Bad randomness. Everything said rests on the independence of trials. A generator with a short or patterned sequence breaks the estimate unnoticeably; which generators will do is in the article on random number generators.

Where it works

• Film and games. Path tracing computes the illumination of each pixel as a mean over random light paths, and the noise on an unfinished frame is Monte Carlo error: halving it takes four times as many rays. This is how, for example, Cycles, Blender’s rendering engine, works. Splitting a pixel into cells with one random ray in each is exactly stratification.
• Finance. Portfolio risk, the price of complex options, payment scenarios: wherever the outcome depends on many random quantities at once.
• Physics. Transport of neutrons and photons, molecules, scattering in detectors: the method began with this problem.
• Geometry and engineering. Volumes and areas of bodies given only by an “inside or outside” test, tolerances of an assembly of many parts.

History

The idea came to Stanisław Ulam in 1946 over a game of patience: the probability that it comes out is easier to estimate by dealing it a hundred times than to compute by combinatorics. Together with John von Neumann he applied it to neutron transport, and their colleague Nicholas Metropolis suggested the name, after the casino. The first paper titled “The Monte Carlo Method” was by Metropolis and Ulam, in 1949. This story is told in detail in the course chapter on randomised algorithms, which also has the method’s best-known experiment: the number π with darts. The interval for a proportion that we used was published by Edwin Wilson in 1927, twenty years before the method had a name.

Problems

In the first you have to give an honest interval from the results of an experiment. In the second Monte Carlo is powerless, and the answer has to be computed exactly.