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

2 Responses to “Buffon versus Bertrand in R”

  1. 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!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.