March 14, 2024

Overview

Monte Carlo simulations are powerful tools that leverage randomness to model various statistical phenomena. In honor of Pi (\(\pi\)) Day, we will use this type of method to estimate the value of \(\pi\).

pi
## [1] 3.141593

Estimating \(\pi\) with Monte Carlo Simulation

Imagine you’re engaged in a friendly game of darts, where a circular dartboard is mounted on a square piece of wood of the same diameter. Your friend claims that they can estimate the number \(\pi\) by randomly throwing darts at the board.

You’re skeptical, but you decide to see how it plays out!

Dart Throwing in Action – 100 Darts!

## [1] "Pi Estimate: 3.12"

Dart Throwing in Action – 10,000 Darts!

## [1] "Pi Estimate: 3.1288"

Dart Throwing in Action – 30,000 Darts!

## [1] "Pi Estimate: 3.1544"

Understanding the Circle-Square Area Relationship

As the darts fly, you notice your friend tallying the darts landing inside vs. outside of the circle (but still within the square). They keep mentioning numbers that seem to gradually get closer and closer to \(\pi\).

To understand what your friend is doing, let’s examine the relationship between the unit circle and the unit square. We can start with the familiar formula for the area of a (unit) circle, \(\pi r^2\). Similarly, we can find that the area of a unit square is \((2r)^2\).

The Ratio of the Area of the Circle to the Area of the Square

Since the area of the circle is \(\pi r^2\) and the area of the square is \((2r)^2\), the ratio of their areas can be calculated as follows:

\[\pi r^2 / (2r)^2 = \pi r^2 / 4 r^2 = \pi/4.\]

pi/4
## [1] 0.7853982

If we take this ratio and multiply it by 4, we get an estimate of \(\pi\)!

Estimating \(\pi\)

The more darts your friend throws, the closer their estimate gets to the actual value of \(\pi\).

Let’s Join the Dart Throwing!

# Import necessary libraries
library(grid)
library(plotrix)

# Define helper functions
is_in_circle <- function(x,y) {
  sqrt((x-0.5)^2 + (y-0.5)^2) <= 0.5
}
estimate_pi <- function() {
  4*length(darts[darts])/num_attempts
}
throw_darts <- function(n) {
  num_attempts <<- n  
  dart_positions <<- matrix(runif(2*num_attempts), ncol=2) 
  darts <<- is_in_circle(dart_positions[,1], dart_positions[,2]) 
}

Let’s Join the Dart Throwing!

plot_darts <- function() {
  plot(c(0, 1), c(0,1), type = "n", asp = 1, 
       main="The Dartboard",
       xlab='Unit Square/Circle',
       ylab='',yaxt='n', frame.plot = FALSE)
  rect(0,0,1,1,border='black')
  draw.circle(0.5,0.5,0.5,border='black')
  points(dart_positions[!darts,1], 
         dart_positions[!darts,2], col='grey', pch=20) 
  points(dart_positions[darts,1], 
         dart_positions[darts,2], col='purple', pch=20)  
}
calculate_pi <- function(n) {
  throw_darts(n)
  estimate_pi()
}

Let’s Throw \(\pi \times 10,000\) Darts!

throw_darts(31416); plot_darts(); estimate_pi()

Let’s Throw \(\pi \times 10,000\) Darts!

## [1] 3.148714

Thank You!