# \pi day!

March 14, 2010
By

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

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.

Tags: , , , , , , ,