\pi day!

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)