Skip to content

Monte Carlo methods

Probability · level 8, lesson 7 of 7

A question to start

Draw a square on the ground, 1 metre on each side, and inside it draw a quarter circle of radius 1 metre with its centre at one corner. Now throw grains of rice at the square from far away, so that they land at random all over it.

Count the grains inside the quarter circle and divide by all the grains in the square. Multiply by 4. You get, roughly, \(\pi = 3.14159\dots\)

Rice can calculate \(\pi\). This is not a trick; it is a method used every day by scientists and engineers, called the Monte Carlo method: using random numbers to compute things that are not random at all. This last lesson shows why it works, how accurate it is, and a famous experiment with needles that also finds \(\pi\).

The idea

Probability as area, run backwards

In level 6 you used geometric probability: a point chosen at random in a region lands in a part of it with probability \(\tfrac{\text{area of part}}{\text{area of region}}\). Monte Carlo turns this round: if you can simulate the random point, the observed fraction estimates the area.

Choose a random point \((U, V)\) in the unit square, with \(U\) and \(V\) independent and uniform on \([0, 1]\) (lesson 2). It lands inside the quarter circle when \(U^2 + V^2 \le 1\). The quarter circle has area \(\tfrac14\pi \cdot 1^2 = \tfrac{\pi}{4}\) and the square has area 1, so

\[P(U^2 + V^2 \le 1) = \frac{\pi}{4}.\]

Repeat \(n\) times and let \(H\) be the number of hits. By the law of large numbers (lesson 5), \(\tfrac Hn \to \tfrac{\pi}{4}\), so

\[\hat\pi = 4 \cdot \frac{H}{n}\]

is an estimate of \(\pi\) that gets better as \(n\) grows. (The hat on \(\hat\pi\) means "estimate of".)

0 1 1

In that picture 200 random points were drawn and 159 landed inside, giving \(\hat\pi = 4 \times \tfrac{159}{200} = 3.18\).

How accurate is it?

Each point is a trial with success chance \(p = \tfrac{\pi}{4} \approx 0.785\). So \(H \sim \operatorname{B}(n, p)\), and the estimate \(4\tfrac Hn\) has standard deviation

\[4\sqrt{\frac{p(1 - p)}{n}} \approx \frac{1.64}{\sqrt n}.\]

One computer run gave:

points \(n\) hits \(H\) estimate \(4H/n\) standard deviation \(1.64/\sqrt n\)
100 83 3.32 0.16
1 000 778 3.112 0.052
10 000 7 801 3.1204 0.016
100 000 78 531 3.14124 0.0052
1 000 000 785 697 3.14279 0.0016

Every estimate is within about 2 standard deviations of the true \(\pi\), just as the normal approximation says it should be about 95% of the time. And notice the price of accuracy: each extra correct decimal place costs 100 times as many points, because the error shrinks like \(\tfrac{1}{\sqrt n}\). A million points give about two or three decimal places. As a way to compute \(\pi\) this is hopeless (formulas give millions of digits); the point is that the same method works for problems no formula can touch.

Monte Carlo integration

The same idea finds any area, and so any integral. To estimate \(I = \int_0^1 g(x)\,dx\), notice that if \(U\) is uniform on \([0, 1]\) then \(E\big(g(U)\big) = \int_0^1 g(x) \cdot 1\,dx = I\) (lesson 2's formula for an expected value). So

\[I \approx \frac{g(U_1) + g(U_2) + \dots + g(U_n)}{n},\]

the average of \(g\) at \(n\) random points, with an error that shrinks like \(\tfrac{\sigma}{\sqrt n}\) by the central limit theorem.

Example. The bell-curve integral \(\int_0^1 e^{-x^2}\,dx\) has no formula in terms of ordinary functions. Averaging \(e^{-u^2}\) over 10 000 random \(u\) gave 0.7460, with a standard deviation of about 0.002. The true value is 0.74682.

The real power appears in many dimensions. A grid with 10 points along each axis needs \(10^{d}\) points in \(d\) dimensions: \(10^{20}\) for 20 dimensions, which is impossible. Monte Carlo's error is \(\tfrac{\sigma}{\sqrt n}\) whatever the dimension. That is why it is used to price financial options, render realistic light in films (tracing random light paths), model particles in physics, and plan the paths of spacecraft through uncertain conditions.

The story goes that the method got its name in the 1940s from the Monte Carlo casino in Monaco, when Stanislaw Ulam, John von Neumann and colleagues used random sampling for calculations in nuclear physics.

Buffon's needle

A much older experiment, posed by the French naturalist Georges-Louis Leclerc, Comte de Buffon, in the 1700s: a floor has parallel lines 1 unit apart. Drop a needle of length 1 at random. What is the chance it crosses a line?

x ½ sin θ θ spacing 1 needle length 1

Describe where the needle lands with two numbers:

  • \(x\), the distance from the needle's centre to the nearest line: uniform on \(\left[0, \tfrac12\right]\);
  • \(\theta\), the angle between the needle and the lines: uniform on \(\left[0, \tfrac{\pi}{2}\right]\) (by symmetry, angles beyond \(\tfrac{\pi}{2}\) behave the same).

The half-needle on the side of that line reaches a height \(\tfrac12\sin\theta\) above (or below) the centre. So the needle crosses exactly when

\[x \le \tfrac12 \sin\theta.\]

The pair \((\theta, x)\) is a random point in a rectangle of width \(\tfrac{\pi}{2}\) and height \(\tfrac12\), area \(\tfrac{\pi}{4}\). The crossing region is the area under the curve \(x = \tfrac12\sin\theta\):

\[\int_0^{\pi/2} \tfrac12 \sin\theta\,d\theta = \tfrac12\Big[-\cos\theta\Big]_0^{\pi/2} = \tfrac12(0 + 1) = \tfrac12.\]

So

\[P(\text{cross}) = \frac{1/2}{\pi/4} = \frac{2}{\pi} \approx 0.6366.\]

\(\pi\) appears because the needle's angle is uniformly spread around a circle. Turning it round, dropping \(n\) needles and counting \(C\) crossings gives the estimate \(\hat\pi = \tfrac{2n}{C}\). A computer dropped 10 000 virtual needles and got 6266 crossings: \(\hat\pi = \tfrac{20\,000}{6266} \approx 3.19\). The crossing fraction 0.6266 is about 2 standard deviations below 0.6366 (the standard deviation is \(\sqrt{0.6366 \times 0.3634 / 10\,000} \approx 0.0048\)): an unluckier run than most, but well within what chance allows. Buffon's needle is a slow way to find \(\pi\), but a beautiful link between geometry and probability.

Pseudo-random numbers

Computers do not toss coins. They use pseudo-random number generators: fixed formulas whose output looks random and passes statistical tests, but is completely determined by a starting value called the seed. That is useful: running again with the same seed repeats the same "random" numbers, so results can be checked. The simulations in this level (the dice average, the random walk, the points above) were all made this way.

Worked examples

1. A simulated probability. To check the gambler's-ruin answer \(\tfrac{3}{10}\) from lesson 6, a computer played 100 000 fair games from 3 dollars with goal 10, and won 29 853. Estimate: 0.2985. The standard deviation of the estimate is \(\sqrt{0.3 \times 0.7 / 100\,000} \approx 0.0014\), so 0.2985 is about 1 standard deviation from 0.3: consistent.

2. How many points for a target accuracy? To estimate \(\pi\) by the quarter circle with a standard deviation of 0.001: \(\tfrac{1.64}{\sqrt n} = 0.001\), so \(\sqrt n = 1640\) and \(n \approx 2.7\) million points.

3. Estimating an area. An odd-shaped lake lies inside a 2 km by 3 km rectangle on a map. A computer picks 5000 random points in the rectangle and finds 1850 in the lake. Estimated area: \(6 \times \tfrac{1850}{5000} = 2.22\) km². Standard deviation: \(6\sqrt{0.37 \times 0.63 / 5000} \approx 0.041\) km².

4. A probability with no easy formula. What is the chance that three random points in a unit square form an obtuse triangle? There is an exact answer, but it is hard; a simulation of a million triangles gives it to about three decimal places in a second. Simulation is often the fastest way to check a difficult calculation, or to find an answer when no one knows how to calculate it.

Try it

  1. In a quarter-circle experiment, 7 900 of 10 000 points land inside. Estimate \(\pi\) and give the standard deviation of the estimate.
  2. Explain why multiplying the number of points by 4 only halves the error.
  3. Describe a Monte Carlo estimate of \(\int_0^2 x^3\,dx\) using uniform random numbers on \([0, 2]\). What is the exact value?
  4. In Buffon's experiment with needle length equal to line spacing, 500 needles give 320 crossings. Estimate \(\pi\).
  5. Check the integral \(\int_0^{\pi/2} \tfrac12\sin\theta\,d\theta = \tfrac12\) by finding the area under the curve with a sketch: is \(\tfrac12\) a reasonable fraction of the rectangle?
  6. Describe how to estimate the chance that the sum of three dice is at least 15 by simulation. Then find the exact answer by counting, and say how many simulated rolls you would need for an error of about 0.005.
  7. A needle of length \(\ell\) less than the spacing 1 is used instead. Show that the crossing chance becomes \(\tfrac{2\ell}{\pi}\).
  8. Why is it useful that a pseudo-random generator gives the same numbers from the same seed?
Answers
  1. \(\hat\pi = 4 \times 0.79 = 3.16\). Standard deviation \(\approx \tfrac{1.64}{\sqrt{10\,000}} = 0.0164\). (The true \(\pi\) is about 1.1 standard deviations away.)
  2. The error is proportional to \(\tfrac{1}{\sqrt n}\); multiplying \(n\) by 4 multiplies \(\sqrt n\) by 2, so the error is halved.
  3. If \(U\) is uniform on \([0, 2]\) (density \(\tfrac12\)), then \(E(U^3) = \int_0^2 x^3 \cdot \tfrac12\,dx\), so the integral is \(2E(U^3)\). Estimate: twice the average of \(U_i^3\) over many random \(U_i\). Exact: \(\tfrac{2^4}{4} = 4\).
  4. \(\hat\pi = \tfrac{2 \times 500}{320} = 3.125\).
  5. The curve rises from 0 to \(\tfrac12\) over the width \(\tfrac{\pi}{2} \approx 1.57\); the rectangle has area about 0.785. The region under the curve is a bit more than half the rectangle, since the sine curve bulges upward. \(\tfrac12\) out of 0.785 is about 64%: reasonable.
  6. Simulate three dice many times and record the fraction with sum at least 15. Exact: sums 15, 16, 17, 18 have 10, 6, 3 and 1 outcomes, total 20 out of 216, so \(p = \tfrac{20}{216} \approx 0.0926\). The standard deviation of the estimate is \(\sqrt{\tfrac{p(1 - p)}{n}}\); for 0.005 we need \(n \approx \tfrac{0.0926 \times 0.9074}{0.005^2} \approx 3361\) rolls.
  7. Now the half-needle reaches \(\tfrac{\ell}{2}\sin\theta\), and the crossing area is \(\int_0^{\pi/2} \tfrac{\ell}{2}\sin\theta\,d\theta = \tfrac{\ell}{2}\). Divided by the rectangle's area \(\tfrac{\pi}{4}\): \(\tfrac{2\ell}{\pi}\). (The condition \(\ell \le 1\) makes sure the needle can cross at most one line.)
  8. Results can be repeated exactly: someone checking your work, or you fixing a bug, can rerun the simulation and get the same numbers. Changing the seed gives a fresh, independent-looking run.

Common mistake

Expecting a simulation to be exact. A Monte Carlo answer is an estimate with a random error of about \(\tfrac{\sigma}{\sqrt n}\). Always report how many trials you used, and ideally the standard deviation, so others know how much to trust it.

Using too few trials. A hundred trials pin a probability down to only about \(\pm 0.05\) to \(\pm 0.1\). Before concluding anything surprising from a simulation, work out its standard deviation.

Where it leads