(This article was first published on

**Xi'an's Og » R**, and kindly contributed to R-bloggers)**F**ollowing my earlier post on Buffon’s needle and Bertrand’s paradox, above are four outcomes corresponding to four different generations (among many) of the needle locations. The upper right-hand side makes a difference in the number of hits (out of 1000). The R code corresponding to this generation was made in the métro, so do not expect subtlety:

#Several ways of throwing a needle at "random" L=0.35 #half-length of the needle D=20 #length of the room N=10^3 numbhits=function(A,B){ sum(abs(trunc(A[,2])-trunc(B[,2]))>0)} par(mfrow=c(2,2),mar=c(1,1,1,1)) #version #1: uniform location of the centre U=runif(N,min=L,max=D-L) #centre O=runif(N,min=0,max=pi) #angle C=cbind(runif(N,0,D),U) A=C+L*cbind(cos(O),sin(O)) B=C-L*cbind(cos(O),sin(O)) plot(C,type="n",axes=F,xlim=c(0,D),ylim=c(0,D)) for (t in 1:N) lines(c(A[t,1],B[t,1]),c(A[t,2],B[t,2]),col="steelblue") for (i in 1:(D-1)) abline(h=i,lty=2,col="sienna") title(main=paste(numbhits(A,B),"hits",sep=" ")) #version #2: uniform location of one endpoint U=runif(N,min=2*L,max=D-(2*L)) #centre O=runif(N,min=0,max=2*pi) #angle A=cbind(runif(N,0,D),U) B=A+2*L*cbind(cos(O),sin(O)) plot(A,type="n",axes=F,xlim=c(0,D),ylim=c(0,D)) for (t in 1:N) lines(c(A[t,1],B[t,1]),c(A[t,2],B[t,2]),col="steelblue") for (i in 1:(D-1)) abline(h=i,lty=2,col="sienna") title(main=paste(numbhits(A,B),"hits",sep=" ")) #version #3: random ray from corner O=runif(N,min=0,max=pi/2) #angle U=L+runif(N)*(D*sqrt(1+apply(cbind(sin(O)^2,cos(O)^2),1,min))-2*L) #centre C=cbind(U*cos(O),U*sin(O)) A=C+L*cbind(cos(O),sin(O)) B=C-L*cbind(cos(O),sin(O)) plot(C,type="n",axes=F,xlim=c(0,D),ylim=c(0,D)) for (t in 1:N) lines(c(A[t,1],B[t,1]),c(A[t,2],B[t,2]),col="steelblue") for (i in 1:(D-1)) abline(h=i,lty=2,col="sienna") title(main=paste(numbhits(A,B),"hits",sep=" ")) #version #4: random ray from corner O=runif(N,min=0,max=pi/2) #angle U=runif(N)*(D*sqrt(1+apply(cbind(sin(O)^2,cos(O)^2),1,min))-2*L) #centre A=cbind(U*cos(O),U*sin(O)) B=A+2*L*cbind(cos(O),sin(O)) plot(A,type="n",axes=F,xlim=c(0,D),ylim=c(0,D)) for (t in 1:N) lines(c(A[t,1],B[t,1]),c(A[t,2],B[t,2]),col="steelblue") for (i in 1:(D-1)) abline(h=i,lty=2,col="sienna") title(main=paste(numbhits(A,B),"hits",sep=" "))

When running the R code for 10⁶ iterations, the approximations to π based on the standard formula are given by

[1] 3.194072 [1] 3.140457 [1] 3.213596 [1] 3.210177

Filed under: R, Statistics Tagged: Bertrand's paradox, Buffon's needle, R, sigma-algebra

To

**leave a comment**for the author, please follow the link and comment on his blog:**Xi'an's Og » 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...