Monte Carlo: estimating by chance

What you will be able to do

After this lesson you can

• estimate a probability by simulation,
• estimate a number, such as $\pi$, from random points, and
• say how the error of such an estimate falls with the number of trials.

Chance as a calculator

The Monte Carlo method, named after the casino, uses random trials to find an answer. The idea is simple: to find the probability of an event, do the experiment many times and count how often the event happens. By the law of large numbers, the share of successes comes near the probability.

This is useful when the calculation is hard. In physics, Monte Carlo methods follow the paths of particles through a detector, and the neutrons in a reactor. In finance, they price options. In every case, the idea is the one you see here.

Two dice

Predict first. Of the 36 pairs of two dice, 21 have the sum 7 or more. What probability do you expect? Then run it.

The simulation says 0.584 and the exact value is $21/36 = 0.583$. For the sum equal to 7: 0.166, and the exact value is $6/36 = 0.167$. The simulation agrees with the exact answer to about three decimals. The simulation is easy to adapt: for three dice, or for a sum above 15, you change one line, while the exact count by hand gets harder.

The birthday problem

Predict first. In a group of 23 people, what is the chance that two have the same birthday? Guess before you run it.

How the program works: each row of days is one group, with a random birthday (0 to 364) for each person. days.sort(axis=1) sorts each row, so that equal birthdays are next to each other, and np.diff(days, axis=1) == 0 finds a neighbour that is equal. .any(axis=1) asks whether a group has at least one such pair, and .mean() is the share of such groups.

Estimating pi

You can find $\pi$ by throwing darts. Imagine a square of side 1, and inside it a quarter of a circle of radius 1. The square has the area 1, and the quarter circle has the area $\pi/4$. If you throw points at random into the square, the share that land inside the quarter circle is $\pi/4$. So $\pi \approx 4 \times$ that share.

The curve is a close formula. The code checks it the Monte Carlo way: make 10,000 random groups and count.

How fast does the error fall?

Predict first. If you use 100 times more points, how much more accurate will the answer be?

A hundred times more points makes the answer only about 10 times more accurate: the printed ratios are 9.8 and 8.1. The error of a Monte Carlo estimate falls as

where $N$ is the number of trials. This is slow. To get one more correct digit, you need 100 times as many trials. So Monte Carlo is not the best way to find $\pi$: there are quick formulas for that. It is the best way to solve problems that have no quick method, such as problems in many dimensions. And the $1/\sqrt{N}$ law is also true for the errors of measurements, as the last lesson shows.

Common mistakes

• Using too few trials. The estimate has a large scatter.
• Trusting one run. Repeat the experiment to see how much the answer changes.
• Using the same seed and thinking that two runs are independent. They are not.
• Forgetting that the share of True is the mean of the comparison.
• Expecting the error to fall quickly. It falls as one over the square root of $N$.

Practice

Start with two quick questions. Each New numbers gives you a fresh one, so you can practise as long as you like.

Then the labs and the problems:

1. Warm-up. Open the lab Estimate a probability. You estimate a probability and $\pi$ with your own simulation.
2. Build. Open the problem Probability of a sum and the problem Estimate pi.

Remember these

• Monte Carlo: estimate a probability as the share of random trials in which the event happens.
• $\pi \approx 4 \times$ the share of random points of a square that fall inside the quarter circle.
• The error falls as $1/\sqrt{N}$: 100 times more trials give 10 times more accuracy.
• Repeat the experiment to find out how large the error is.

Check yourself

Try each one on your own first. If you are stuck, open a hint. Then check with Show the answer.

1. What is the probability that a die shows an even number, and how can you estimate it?

2. In the pi experiment, 786 of 1000 points are inside. What is the estimate of $\pi$?

3. By how much do you reduce the error if you use 10 000 times more trials?

4. Why is the result for 23 people in the birthday problem so large?