How to make pi ... using R

There's another irresistible post over at Dot Physics, this time on a nifty way to estimate the value of pi using random numbers. Check out that post for details, then hop back over here!

(brief pause...)

Ok, you're back! And you loved the post - awesome.  But didn't it feel like it was lacking a little... pie?  I too had the very same feeling, so I wrote my own version of Rhett's simulations (using R) but with some fanciful graphics just to jazz things up a bit.

Figure 1: The fraction of points that fall "on the pie" (left) give a reasonable 

approximation of π=3.14159... Lines on the right show π and our estimate.

If you have R installed on your computer, run the code below for the animated version.

## Estimating pi... with pie.
##

## Initialize number of points to use (5 reps of each)
Ns=sort(runif(50, min=1, max=10000));   
pi.est=Ns*NA;  ## place holder for estimates

## Open up a window with room for two figures.
X11(14,7); par(mfrow=c(1,2))


## Loop through the different numbers of points
for(i in 1:length(Ns)) {
   N=Ns[i];  ## Current number of points
   XY=matrix(runif(2*N),N,2) ## randomly place each point.
   pie  = sqrt(XY[,1]^2+XY[,2]^2)<1; ## points on the pie
   plate=!pie;                       ## points off the pie

 ## Draw the pie...
   plot(c(0,1),c(0,1),type="n", xlab="", ylab="")
   xs=seq(0,1,length=100)
   polygon(c(xs,rev(xs)), c(sqrt(1-xs^2),0*xs), 
           col="khaki", border=NA)

 ## Draw the points on the pie, and outside the pie
   symbols(XY[pie,], circles=rep(0.01,sum(pie)), 
        fg="khaki1", bg="khaki1",
        inches=FALSE, add=TRUE,)
   points(XY[plate,], col="lightgray", pch=20, cex=2)

 ## It's pie, so why not slice it?  
   points(c(0,sqrt(1/2)), c(0,sqrt(1/2)), type="l")

 ## Fruit pie? Lets decorate...
   points(c(0.6, 0.7)*cos(3*pi/8+.06),  ## using polar coords
          c(0.6, 0.7)*sin(3*pi/8+.06),  ## to make 4 vents.
          type="l", lwd=6, col="red3")
   points(c(0.6, 0.7)*cos(3*pi/8-.06), 
          c(0.6, 0.7)*sin(3*pi/8-.06), 
          type="l", lwd=6, col="red3")
 
   points(c(0.6, 0.7)*cos(1*pi/8+.06), 
          c(0.6, 0.7)*sin(1*pi/8+.06), 
          type="l", lwd=6, col="red3")
   points(c(0.6, 0.7)*cos(1*pi/8-.06), 
          c(0.6, 0.7)*sin(1*pi/8-.06), 
          type="l", lwd=6, col="red3")

 ## Plot estimates of pi in a second window... 
   pi.est[i]=4*sum(pie)/N; 
   title(bquote(hat(pi) == .(pi.est[i])))  ## See '?plotmath'
   plot(Ns, pi.est, pch=20, xlab="N", ylab=bquote(hat(pi)), 
        xlim=range(Ns), ylim=c(0,4), color="darkred", 
        main=expression(paste("Estimates of ", pi)))
   Sys.sleep(0.25); ## pause briefly between pies.
} # end for loop


## Evaluate our estimates in the "large N" range
y=subset(pi.est, Ns>median(Ns));
abline(h=pi);  ## draw lines at y=pi, y=mean(pi.est)
abline(h=mean(y), lty=2, col="red");

text(median(Ns), 2, paste("N>500 :: mean =",signif(mean(y),5), 
                          ", sd =",signif(sd(y),3)))

Oh, and before anyone else pounces on the obvious, you can of course make even better pi using R via...

> pi
[1] 3.141593
>