Estimating Pi Using Monte Carlo Simulations

Pi (3.141593) is one of the few magical numbers in Mathematics that we often trust, accept and use in our calculations. However, you might be curious to know where it comes from. Pi can be obtained analytically which gives us a value equal to 3.141593 but here we're going to find the value of pi numerically by running Monte Carlo simulations.

As you might remember from primary or elementary school, the formula to obtain the area of a circle is:

A = \pi \times r^2

where $r$ is the radius of the circle. So if we have a circle with $r=1$ , then its area is just $\pi$ !

Now imagine we have a circle inside a square with sides equal to 2 and if we look at a corner of it, we see that $\frac{1}{4}$ of the area of the circle is inside a square with sides equal to 1. So we know the area of that slice of circle is $\frac{\pi}{4}$ .

So if we denote the area of the quarter of the circle as $Q$ and we know:

\pi = 4 \times \frac{Q}{r^2}

Now, using Monte Carlo simulations, we generate a large number, lets say a million, $x$ and $y$ coordinates that are uniformly distributed $Unif(0, 1)$ and using the pythagorus formula we find their distances from the centre of the circle

r' = \sqrt{x^2 + y^2}

If our ${r'}$ is smaller than the radius of the circle, $r=1$ , then it's inside the circle and vice versa. If we run enough simulations, then the proportion of points that are inside the circle to the total number of trials would essentially give us the area of that slice of the circle, $Q$ .

For a circle with $r=1$ , if we multiply that area by 4, using the $\pi$ formula above, we get our estimate for $\pi$ .

Plot of 500,000 simulations

Plot of 1,000,000 simulations

Code


library(ggplot2)
trials <- 1000000

r <- 1
x <- runif(trials, 0, r)
y <- runif(trials, 0, r)

distance_from_center <- sqrt(x^2 + y^2)
inbounds <- distance_from_center < r


ggplot(data.frame(x, y, inbounds), aes(x, y, color=inbounds)) +
  theme_bw() +
  ggtitle(paste(round(trials), 'Trials')) +
  guides(color=FALSE) +
  geom_point(size=0.002) +
  ggtitle(paste(trials, 'Trials'))

pi_estimate <- 4 * sum(inbounds)/trials
pi_estimate

error = 1 - pi_estimate/pi
error