Estimating Pi Using Monte Carlo Simulation in R

Everyone knows that \(\pi\) is something about 3.14…

Lets say that we want to estimate \(\pi\) but do not know any formula to do that. This is a problem! Fortunately, we know something about areas of square and circle:

Area of square \[ A = a^2 \]

Area of circle \[ A = \pi r^2 \]

Suppose that we have a circle that has a radius that is half the length of the side of the square. Now we can write the area of the square as:

\[ A = (2r)^2 \]

\[ A = 4r^2 \]

Now, if we know both areas we can estimate \(\pi\) using the formula:

\[ \frac{A_c}{A_s} = \frac{\pi r^2} {4r^2} = \frac{\pi}{4} \]

Ok, but how can we get both areas without using any mathematical formulas? The answer is simple: we just need to use Monte Carlo method and approximate both areas using some random process.

Let’s draw a sample of 1000 random points that lie inside the 2 by 2 square:

x <- runif(1000, -1, 1) # 1000 random numbers between -1 and 1
y <- runif(1000, -1, 1)

And plot them:

plot(x,y, pch = 20)

Some of them lie inside a circle with radius one. The question is which? If we know that \(r=1\) we can compute distance from the origin \((0, 0)\) using the formula:

\[ d = \sqrt{x^2 + y^2} \]

If it is smaller than radius (\(1\)) it means that the point lies inside the circle!

colors = c() # create vector with colors
for (i in 1:length(x)){
  dist <- sqrt(x[i]**2 + y[i]**2) # formula for distance from the origin
  if (dist < 1){
    color <- "red" # red if the point lies inside the circle 
  } else {
    color <- "black" # black if it does not
  }
  colors <- c(colors, color)
}

plot(x,y, pch = 20, col = colors)

Now, to approximate \(\frac{\pi}{4}\) we need to compute the ratio between red points and all points!

count_red <- table(colors)["red"]
count_red / 1000

##   red 
## 0.789

To get \(\pi\) we just need to multiply it by four (because the ratio is \(\frac{\pi}{4}\):

(count_red / 1000) * 4

##   red 
## 3.156

Pretty bad, huh? Maybe we have too few points. Lets try with 100 000:

n_sim <- 100000
x <- runif(n_sim, -1, 1)
y <- runif(n_sim, -1, 1)

in_circle <- 0 # we will use simple counter

for (i in 1:length(x)){
  dist <- sqrt(x[i]**2 + y[i]**2)
  if (dist < 1){
    in_circle <- in_circle + 1
  } 
}

(in_circle / n_sim) * 4

## [1] 3.14092

This is much better! The basic idea is that by selecting random points from uniform distribution we can approximate ratio between area of the circle and area of the square and by doing that we can estimate \(\pi\)!

To be frank this example has been done to death. Let’s do something more interesting. Instead of 2D shapes we can use 3D shapes!

Volume of cube \[ V = a^3 = (2r)^3 = 8r^3 \] Volume of sphere \[ V = 4/3 \pi r^3 \]

Ratio \[ \frac{V_s}{V_c} = \frac{4/3 \pi r^3}{8r^3} = \frac{4/3\pi}{8} = \frac{\pi}{6} \]

n_sim <- 10000
# Now we have 3 coordinates
x <- runif(n_sim, -1, 1) 
y <- runif(n_sim, -1, 1)
z <- runif(n_sim, -1, 1)
library(rgl) # library to plot in 3D
in_sphere <- 0
colors = c()
for (i in 1:length(x)){
  dist <- sqrt(x[i]**2 + y[i]**2 + z[i]**2) # distance from the origin
  if (dist < 1){
    color <- "red"
    in_sphere <- in_sphere + 1
  } else {
    color <- "black" 
  }
  colors <- c(colors, color)
}
plot3d(x,y,z, pch = 20, col = colors)
(in_sphere / n_sim) * 6

## [1] 3.1746

Enjoy the Pi Day!

Today is Pi Day!

Estimating Pi Using Monte Carlo Simulation in R