The Monte Carlo Pi

In this post will you see, a tangled explanation of the Monte Carlo simulation procedure, finding the value of pi using a geometric concept where area and probability have a relation. Then follows the concept of computing the standard error to estimate the variability created by errors in measurement (or maybe that will be another post).

Pi is the ratio of the circumference of a circle to its diameter. To estimate the value of Pi, there are many numerical estimators, for instance the one you learned in school was 22/7. Here we will use a geometric method to compute the area of a figure using random numbers thus establishing the relation between the area of a figure and the probability of an event happening in that figure as a function of its area.

So the square is of side 2r where r is the radius of the circle inscribed in it. The area of the square is thus 4*r^2 and the area of the circle obviously is pi*r^2.

The ratio of the area of the circle to the area of the square gives us 1/4th the value of pi. That is (pi*r^2)/(4*r^2) = pi/4.

The equation of a circle is (x-h)^2 + (y-k)^2 = r^2 where (h,k) is the co-ordinate of the center point of the circle and r is the radius. For a unit circle (circle of radius 1) at origin; (h,k) is (0,0) and r^2 =1. So the equation of this circle is x^2+y^2=1. That in simple english means that every point (x,y) co-ordinate that satisfies the relation x^2+y^2=1 lies along the circumference of this unit circle.

So if I were to randomly throw darts at the square, that could be simulated by generating uniformly distributed random points (x,y) co-ordinates that lie between -1 and 1. Since the centre of the circle is 0,0 and the circle is unit radius, the side of the square is 2 units in x-co-ordinate (-1 to 1) and the same for the y-co-ordinate. So if I generate x,y co-ordinates with uniform random numbers ranging from -1 to 1 they will all lie within the square.

The number of darts that would land in the circle is proportional to its area. If it were a very small circle, much lesser darts would land in the circle and if it circumscribed the square, all the darts that landed in the square would land in the circle. So you get a fair idea of what I mean when I say that the area of the circle gives a measure of the probability that a dart would land within the circle.

Now the ratio of the number of darts in the circle divided by the number of darts in the square will give us 1/4th the value of pi. However, how do we estimate the number of points that are within the unit circle? All the randomly generated (x,y) pairs that satisfy the condition x^2+y^2 <=1 are within the unit circle.

If you don't believe me, here is some R code that will do it for you.

x<-runif(100000,min=-1,1)
y<-runif(100000,min=-1,1)
z=x^2+y^2 

plot(x[which(z>1)],y[which(z>1)],pch=".",col="magenta")

points(x[which(z<=1)],y[which(z<=1)],col="cyan",pch=".")
 

First let's see where all the randomly generated (x,y) co-ordinates with a value greater than 1 for the function x^2+y^2 go in the plot. As you can clearly see, all the points in magenta have fallen outside the circle, this means that if you plot the other points in cyan they will all be inside the circle. So our uniformly distributed points have more or less fallen into shapes proportional to their respective areas. Now if we take the ratio of the number of points to have fallen within the circle and within the square (that is the total number of points generated) and multiply this number by 4 we will get an approximation of the value of pi.

So this gets you motivated for the next part which is simple and where we obtain an estimate of pi using simply the ratio of the counts of these points. This value of pi is close but it isn't exactly pi (3.1415926), which means that our measure of pi has some uncertainty, so we go about trying to quantify this uncertainty, since what we did was a stochastic simulation we need to perform more runs of the same simulation so that we get an estimate of how much does the value of pi oscillate around the real value.

length(which(z<=1.0000000))/length(z))*4.0000000

[1] 3.14112

Now we repeat this process say 1000 times and we store the value of pi obtained for each iteration in a vector. After the 1000 runs are over, we look at the distribution of the value of pi over the course of 1000 runs. So now the error in measurement is the standard deviation of the distribution of pi. You can see in the histogram styled density plot that the mean value of the density which is the peak, coincides with the point on the x-axis where the real value of pi must be. This concept is similar to how you estimate a statistic like the Standard Error of Mean (SEM), so you repeatedly take values of the statistic by random sub-sampling, then the distribution of the sub-sampled statistics gives you the standard deviation of the statistic which if it turns out to be the mean is called SEM.

pi<-vector()

pb<-txtProgressBar(1,1000,style=3)
for(i in 1:1000)
{
setTxtProgressBar(pb,value=i)
x<-runif(1000000,min=-1,1)
y<-runif(1000000,min=-1,1)
z=x^2+y^2
pi[i]=(length(which(z<=1.0000000))/length(z))*4.0000000
}

# The mean value of this distribution of "pi-es" will be the value that is closest to pi

 plot(density(pi),col="red",type="h")
 
mean(pi)
3.141579 # which is pretty close to the real value
 

> sd(pi)
[1] 0.001547173

# is the measure of our error in the sampling statistic that we used to estimate the value of pi, we can improve the estimate of the value of pi by increasing both the number of points used as darts in the simulation and also the number of iterations for which we repeat the complete process.