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 is the experiment tha Buffon performed using a needle of length,. 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)

*Related*

To

**leave a comment** for the author, please follow the link and comment on his blog:

** Stats raving mad » R**.

R-bloggers.com offers

**daily e-mail updates** about

R news and

tutorials on topics such as: visualization (

ggplot2,

Boxplots,

maps,

animation), programming (

RStudio,

Sweave,

LaTeX,

SQL,

Eclipse,

git,

hadoop,

Web Scraping) statistics (

regression,

PCA,

time series,

trading) and more...

If you got this far, why not

__subscribe for updates__ from the site? Choose your flavor:

e-mail,

twitter,

RSS, or

facebook...

**Tags:** Buffon's Needle, code, day, estimate, pi, Probability, R, statistics