Estimating Pi using Monte Carlo simulations

Published on 27 Aug 13
monte-carlo

A Monte Carlo simulation is a method of finding an answer to a problem by repeatedly and sampling something at random. They are primarily suited to problems where a deterministic algorithm might not be practical or even feasible.

From a financial perspective, Monte Carlo simulations are an effective way of modelling behaviours that have a degree of uncertainty, which can range from calculating option prices to running risk simulations.

How to find Pi

If we consider a unit circle (a circle with a radius of 1) inscribed within a square (making the square 2 units wide) we can show that the ratio of the area of the circle to that of the square can be used to find the value of Pi.

The area of a square with sides of length \(x\) is obviously \(x*x\) , and we know the area of a circle with a radius \(r\) is:

\begin{equation*} A_c = \pi r^2 \end{equation*}

We can create the following ratio between the two areas, and substitute our values in.

\begin{equation*} R = \frac{\pi r^2}{x^2} = \frac{\pi}{4} \end{equation*}

So now we can see that by finding a value for \(R\) we can find a value for Pi.

\begin{equation*} \pi = 4 R \end{equation*}

Using Monte Carlo simulations can provide a simple way of calculating a value for Pi.

Constructing a simulation

Monte Carlo simulations rely on random sampling, and in this simulation we're going to randomly pick a coordinate in the original square. By keeping track of how many times we find a point that falls inside the circle, we can start to approximate \(R\) .

The image above illustrates what is happening in this simulation. The ratio of the number of red dots to blue dots approximates the ratio of area of the circle to the area of the square, which gives us \(R\) .

from random import random
from math import pow, sqrt

simulations = 5000
hits = 0.0

for x in xrange(simulations):
    # random coordinates
    x = random()
    y = random()
    # check distance from (0,0)
    distance = sqrt(pow(x,2) + pow(y,2))

    # if distance < 1 then it's inside the circle
    if distance <= 1.0:
        hits = hits + 1.0

# Now we have an approximation for R: (hits/simulations)
pi = (hits/simulations) * 4.0

print pi

Results

Simulations	Pi
50	3.2000
250	3.0240
1000	3.1960
10000	3.1568
50000	3.1487

As you can see from the table above, our estimations of Pi eventually start to converge on the real value of Pi. However, increasing the number of simulations increases the time it takes to complete.