Buffon versus Bertrand in R
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
Xi'an's Og 
August 31, 2014 at 6:21 am
Hi,
I know this is like 3 years late, but I thought I should post anyway. I’m putting together a poster on Monte Carlo methods for high school students to accompany an activity on Buffon’s needle, and I really liked your 4 methods of approach, and I’ll be adapting them for a section of the poster.
However, I do think you’re making some systematic errors. For example, in the first one, by restricting the centre to be between l and d-l, you’re completely discounting needles that are within the board (that can be 0.01 away from the bottom but have an angle such that they don’t cross, say 0 degrees). As a result, you underestimate the percentage of crosses
Rewriting your code using floor instead of trunc (for negatively placed needles) in numbhits (can be interpreted as rethrowing them if center doesn’t lie in the range etc) and a proper uniform distribution 0,d, for 10^6 trials, you do get quite close to pi. Similarly for part 2.
U=runif(N,min=0,max=D) #centre
Similarly for 3 and 4, I was a bit confused where the range equation came from for the centres. For each angle (on a square board), it’s range is given by max(d/sin(O), d/cos(O)) obtaining instead:
U=runif(N)*(D*apply(cbind(1/sin(O),1/cos(O)),1,max)) #centre
#Alternatively mirror in pi/4 and stick with one.
Again, running this code now, (without the obvious circle range) now gives close to pi when run 10^6 times ~ 3.144261
Overall, I don’t want to appear to be too critical (and I could be wrong), but I’m quite sure with my results.
Either way, thanks, the visualisations were superb and the code was easy to read and modify, and I’m certainly going to discuss the angle origin variant (3 and 4) and maybe introduce some white noise on it’s final position!
August 31, 2014 at 9:59 am
Thanks for the comments and the corrections for the boundary effects! No debate about this.