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

Following 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