Site icon R-bloggers

Buffon versus Bertrand in R

[This article was first published on Xi'an's Og » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
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

To leave a comment for the author, please follow the link and comment on their blog: Xi'an's Og » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.