Approximating Pi via the Monte Carlo Method

Background

This is the typical introductory application of the Monte Carlo method, and as such you can find simlar examples all over the place. That said I’m learning to build Shiny apps and thought this would be a neat little first attempt. I’m not writiing up the Shiny app here, just the problem itself. The idea is to approzimate pi by generating uniformly distributed pairs of random (real) numbers on [-1,1], such that they all fall in a square of area = 2 centered at the origin. We then look at the proportion of points landing in the unit circle, and use that to approximate pi. As the number of points goes to infinity, the ratio tends to pi.

Setting up the problem

We know the area of the square is 2. The area of the unit circle is pir^2, and since r=1 this equals pi. Therefore the ratio of the areas is pi/4. That means if the ratio of points falling in the unit circle is (points in circle)/(total points), we can quickly solve for pi such that pi=4*(points in circle)/(total points)

Let’s start with plotting a unit circle You can imagine randomly trowing darts at this plot, and figuring out the ratio of darts landing in vs out of the circle (as long as you only count darts landing in the square bounded by -1,1).

1k points

Instead of throwing darts, we’ll generate uniformly distributed random numbers. Let’s start with 1000 and see what happens.

number_points<-1000

x<-runif(number_points,-1,1)
y<-runif(number_points,-1,1)

points<-data.frame(cbind(x,y))

Let’s plot these points on our previous graph

To check if a point lands in unit circle, we need to know if x^2 + y^2 < 1. So let’s add a column to our data frame with this calculation.

points$unit_circle<-points$x^2 + points$y^2

We can then get the length of the vector were this value is < 1 to know how many points landed in the circle, and compare this to the overall points to get the ratio. Then, we can use our previous formula to approximate pi.

4*length(points$unit_circle[points$unit_circle<1]) / number_points
## [1] 3.18

Not great. So let’s increase the number of points and try a couple more.

10k points

4*length(points$unit_circle[points$unit_circle<1]) / number_points
## [1] 3.1448

100k points

4*length(points$unit_circle[points$unit_circle<1]) / number_points
## [1] 3.14236

5M points

4*length(points$unit_circle[points$unit_circle<1]) / number_points
## [1] 3.141658

Getting better! I’ll leave it here, but you can try for yourself with a greater number of points to see that the approximation will continue to improve.