The famous mathematical constant, Pi - we all remember it from school, mostly as a way to find the circumference or the area of a circle. In most cases, for simplicity, the value of Pi was often rounded to 3.14.
However, Pi is in fact what mathematicians call an irrational number, meaning that it is a number that can't be written as a simple fraction, or in other words it has an infinite, non-repeating number of decimal digits. Move over 3.14 - here is the value of Pi to 100 digits:
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
If you think that's too much to retain, in October 2015 Suresh Kumar Sharma of India recited 70,030 digits of Pi in just over 17h!
The fact that Pi is infinite hasn't stopped mathematicians trying to calculate it to as many digits as possible. In November 2016, Peter Trueb used some fairly chunky computing power to calculate it to 22,459,157,718,361 decimal places.
In this post, we look to approximate Pi by throwing darts at a dart board using Monte Carlo simulation. For this, all we'll need R, some random numbers, a square, and a circle.
The circle and the square give us our dartboard and back-board respectively, with the circle fitting perfectly inside the square with a radius of 1, like this:
- The circle, with radius 1 has an area of Pi
- The square with sides of length 2 has an area of 4
- Therefore, the ratio of the area of the circle, to the area of the square is Pi/4.
If we repeatedly throw darts at the board and they are to randomly land within the bounds of the square, some will land on the square and some will land on the circular dartboard. If these throws are truly random, then the number of darts that land on the dartboard, divided by the number that we threw in total, will be in the ratio described above (Pi/4). If we multiply that number by 4, we have our estimate of Pi!
One last thing before we run some code - how will we know if the dart lands within the circle or not? Well, think of each dart throw having an x co-ordinate and a y co-ordinate, both between -1 and 1 (0,0 would be the centre of the dart board). Using Pythagoras's theorem we can take the square-root of (x co-ordinate squared + y co-ordinate squared) - and if this value is less than 1 then we know the dart landed inside the circle.
Let's do this in R! We'll use a For Loop to run through several iterations, each using a greater number of dart throws to see if we can improve our estimation of Pi.
# create empty data frame to append estimations to pi_estimation_final <- data.frame(darts_thrown=numeric(0),pi_estimate=numeric(0))
#loop through varying numbers of darts, in this case we'll start with 2^2 darts (4) and keep doubling this until 2^25 darts (33.5 million darts)
for (i in 2^(2:25)){ throws = i # draw random numbers for the x and y co-ordinates of the darts x_coordinate <- runif(throws,0,1) y_coordinate <- runif(throws,0,1)
# use the pythagorean theorem dart_landing_point <- sqrt((x_coordinate^2) + (y_coordinate^2))
# count the number of darts that landed within the circle darts_inside_circle <- sum(dart_landing_point < 1)
# estimate pi for this iteration pi_estimation <- data.frame(darts_thrown = throws, pi_estimation = (darts_inside_circle/throws)*4)
# append this iteration to our final table
pi_estimation_final <- rbind(pi_estimation_final,pi_estimation)
}
Let's plot our results to see if we get closer to Pi as we use more darts - we'll add a red line showing what we're aiming for
#plot our estimations plot(pi_estimation_final$pi_estimation, type='b', main = 'Estimating Pi using Monte Carlo simulated dart throwing', xlab = 'Number of Darts Thrown (2^2 through 2^25)', ylab = 'Pi estimation', lwd = 2, pch=21, bg="black") abline(h=3.1415926, col = 'red', lwd = 2)
When we ran this - the iteration that used 2^25 darts estimated Pi as 3.141493 so still out by around 0.0001.
Have fun! If you found this interesting, please do share!