How to make pi ... using R

Thursday, March 11, 2010 at 12:04 PM Bookmark and Share
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
>

7 comments:

Posted by: Anonymous | 12/04/2010 3:08 PM

hi 1st comment sweet

Posted by: matt | 12/04/2010 3:10 PM

but how do you make it like 58/99

Posted by: matt | 12/04/2010 3:23 PM

i can get close 29/9.2

Posted by: Paul | 12/04/2010 3:25 PM

Ah, but can you get strawberry-rhubarb? ;)

Posted by: matt | 12/04/2010 3:27 PM

f***ing closer 29/9.231

Posted by: matt | 12/04/2010 3:35 PM

btw the first comment was me

Posted by: matt | 12/04/2010 3:38 PM

is 3.14159 rounded from 3.141588...

Post a Comment