Estimating π using Monte Carlo simulations

The Monte Carlo method

Monte Carlo methods are a group of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle. They are often used in physical and mathematical problems and are most useful when it is difficult or impossible to use other approaches. Monte Carlo methods are mainly used in optimization problems, numerical integration, and to generate draws from a probability distribution.1

Estimating the value of \(\pi\)

Circle diagram

Consider the circle in the diagram. If we generate a sufficiently large number of \((x, y)\) points in the domain \(-1 \leqslant x \leqslant 1\) and \(-1 \leqslant y \leqslant 1\), then the fraction of points that fall inside the circle will be equal to the ratio between the area of the circle and area of the square.

$$\frac{\textrm{number of points that fall inside the circle}}{\textrm{total number of points}} = \frac{\textrm{circle area}}{\textrm{square area}} = $$
$$ = \frac{\pi R^2}{(2R)^2} = \frac{\pi}{4}$$

A point falls inside the circle if it complies to the following condition:

$$x^2 + y^2 \leqslant R^2 $$

By counting the number of random points that fall within the circle we can estimate \(\pi\) as:

$$\pi = 4 \times \frac{\textrm{number of points that fall inside the circle}}{\textrm{total number of points}}$$

The algorithm:

We can thus devise an algorithm that will:

  1. Initialize a counter of the points that fall inside the circle (circle_points).
  2. Generate a random point \((x, y)\).
  3. Check if the point falls within the circle by evaluating \(x^2 + y^2 \leqslant R^2\).
  4. If the abose condition is true then increment the counter circle_points.
  5. Calculate \(\pi\) .
  6. Repeat from step 2 for the desired number of runs.

In principle, the higher the number of runs the better the estimation will be (closer to the real value of \(\pi\)).

The code

Putting the above algorithm to pratice we get the following code2:

import random

def run_simulation(runs=1000):
    simulation = []
    circle_points = 0
    for i in range(runs):
        # Create a (x, y) point at random
        x, y = random.uniform(-1, 1), random.uniform(-1, 1)

        # Check if the point falls within the circle
        if x**2 + y**2 <= 1:
            circle_points += 1

        # Calculate current pi value
        pi_estimate = 4 * circle_points/(i+1)

        # Save simulation step 
        simulation.append([x, y, pi_estimate])

    # Print final pi value
    print('Result: Pi =', simulation[-1][2])
    return simulation

sim = run_simulation(runs=10000)

Here is the result of the simulation from 1000 runs to 100000 runs. At 100000 runs the estimation (\(\pi = 3.142720\)) is already very close the actual value (\(\pi = 3.141593\)).

Simulation

And if you are curious, here is how the estimation approaches the real value from the first run to 10000th run.

Evolution

A final remark: We do not need to perform the simulation above using a full circle. By selecting a quarter of a circle and the corresponding square we could achieve the same result with a fourth of the computational effort.

Notes

  1. Monte Carlo method, Wikipedia.
  2. Full source code