# \pi day!

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**Stats raving mad » R**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

It’s π-day today so we gonna have a little fun today with Buffon’s needle and of course R. A well known approximation to the value of $latex \pi$ is the experiment tha Buffon performed using a needle of length,$latex l$. What I do in the next is only to copy from the following file the function estPi and to use an ergodic sample plot… Lame,huh?

estPi<- function(n, l=1, t=2) { m <- 0 for (i in 1:n) { x <- runif(1) theta <- runif(1, min=0, max=pi/2) if (x < l/2 * sin(theta)) { m <- m +1 } } return(2*l*n/(t*m)) }

So, an estimate would be…

estPi(2000,l=1,t=2) # 3.267974

Ok, not that great but for the whole scene it’s remarkable good! Now, we set some increasing sample sizes to account for the estimation.

n=8000 r=15 mat=rep(NA,r) size=rep(NA,r) for (i in 1:r) { size[i]<-n*i mat[i]<-estPi(n*i,l=1,t=2) } matrix<-expand.grid(size) matrix[,2]<-mat names(matrix)<-list("n","pi") matrix # n pi #1 8000 3.182180 #2 16000 3.165809 #3 24000 3.135615 #4 32000 3.145581 #5 40000 3.138486 #6 48000 3.144860 #7 56000 3.162412 #8 64000 3.111932 #9 72000 3.097574 #10 80000 3.155072 #11 88000 3.157404 #12 96000 3.144139 #13 104000 3.126597 #14 112000 3.150226 #15 120000 3.136599

Which is the best estimate?

matrix[which.min(abs(matrix[,2]-pi)),] # n pi # 12 96000 3.144139 plot(matrix,type="b");abline(h=pi,col="red",lty=2)

source : [Chiara Sabatti , pdf]

Take a look @

+ Wiki

+ An introduction to geometrical probability: distributional aspects with applications (A. M. Mathai)

To

**leave a comment**for the author, please follow the link and comment on their blog:**Stats raving mad » R**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.