In honour of π day (03.14 – can’t wait until 2015~) , I thought I’d share this little script I wrote a while back for an introductory lesson I gave on using Monte Carlo methods for integration.

The concept is simple – we can estimate the area of an object which is inside another object of known area by drawing many points at random in the larger area and counting how many of those land inside the smaller one. The ratio of this count to the total number of points drawn will approximate the ratio of the areas as the number of points grows large.

If we do this with a unit circle inside of a unit square, we can re-arrange our area estimate to yield an estimate of π!

This R script lets us see this Monte Carlo routine in action:

############################################## ### Monte Carlo Simulation estimation of pi ## ## Author: Corey Chivers ## ############################################## rm(list=ls()) options(digits=4) ## initialize ## N=500 # Number of MC points points <- data.frame(x=numeric(N),y=numeric(N)) pi_est <- numeric(N) inner <-0 outer <-0 ## BUILD Circle ## circle <- data.frame(x=1:360,y=1:360) for(i in 1:360) { circle$x[i] <-0.5+cos(i/180*pi)*0.5 circle$y[i] <-0.5+sin(i/180*pi)*0.5 } ## SIMULATE ## pdf('MCpiT.pdf') layout(matrix(c(2,3,1,1), 2, 2, byrow = TRUE)) for(i in 1:N) { # Draw a new point at random points$x[i] <-runif(1) points$y[i] <-runif(1) # Check if the point is inside # the circle if( (points$x[i]-0.5)^2 + (points$y[i]-0.5)^2 > 0.25 ) { outer=outer+1 }else { inner=inner+1 } current_pi<-(inner/(outer+inner))/(0.25) pi_est[i]= current_pi print(current_pi) par(mar = c(5, 4, 4, 2),pty='m') plot(pi_est[1:i],type='l', main=i,col="blue",ylim=c(0,5), lwd=2,xlab="# of points drawn",ylab="estimate") # Draw true pi for reference abline(pi,0,col="red",lwd=2) par(mar = c(1, 4, 4, 1),pty='s') plot(points$x[1:i],points$y[1:i], col="red", main=c('Estimate of pi: ',formatC(current_pi, digits=4, format="g", flag="#")), cex=0.5,pch=19,ylab='',xlab='',xlim=c(0,1),ylim=c(0,1)) lines(circle$x,circle$y,lw=4,col="blue") frame() #blank } dev.off() ############################################## ##############################################

The resulting plot (multi-page pdf) lets us watch the estimate of π converge toward the true value.

At 500 sample points, I got an estimate of 3.122 – not super great. If you want to give your computer a workout, you can ramp up the number of iterations (N) and see how close your estimate can get. It should be noted that this is not an efficient way of estimating π, but rather a nice and simple example of how Monte Carlo can be used for integration.

In the lesson, before showing the simulation, I started by having students pair up and manually draw points, plot them, and calculate their own estimate.

If you use this in your classroom, drop me a note and let me know how it went!

I wrote some code to do this too, thought I’d share it incase you are interested. There are no loops or data.frames in this code so it runs quite a lot faster;

n <- 5e4L

x <- runif(n, -1, 1)

y <- runif(n, -1, 1)

isInCircle <- x^2 + y^2 <= 1

est <- 4 * sum(isInCircle) / n

plot(x, y, type = "n")

points(x[isInCircle], y[isInCircle], col = "red", pch = ".")

points(x[!isInCircle], y[!isInCircle], col = "black", pch = ".")

title(substitute(paste("estimate of ", pi, " = ", x), list(x = est)))

Thanks, Jeff!

I wrote it with loops for illustrative purposes since the script plots out a multipage pdf which you can flip through and see the result of each additional point.

I also noticed on r-bloggers that someone seems to invariably do a post on this every year on pi day. Hope you had a good one!

Pingback: Pi in Pubmed « Jean-Etienne's blog

Pingback: Pi day | God plays dice

Reblogged this on lava kafle kathmandu nepal.