*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.

Oh, and before anyone else pounces on the obvious, you can of course make even better## 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)))

*pi*using R via...

> pi [1] 3.141593 >

## 7 comments:

hi 1st comment sweet

but how do you make it like 58/99

i can get close 29/9.2

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

f***ing closer 29/9.231

btw the first comment was me

is 3.14159 rounded from 3.141588...

## Post a Comment