The Buffon’s Needle Problem

Attached is a short write-up on the very interesting geometric probability problem commonly referred to as the Buffon’s Needle problem.  The solution gives a general outline for a Monte Carlo method of approximating pi. neato!Buffons Needle\

An R implementation of the monte carlo simulation is:

#Enter the desired number of trials as n

Buffon_Needle<-function(n){
	L=2;D=L;

	thetas<-c(rep(0,n))
	for(i in 1:n){thetas[i]=runif(1,0,pi/2)}

	X<-c(rep(0,n))
	for(i in 1:n){X[i]=runif(1,0,L/2)	}

	COS<-c(rep(0,n)); COS=cos(thetas)
	Compare<-c(rep(0,n))
	Compare=(X/(2/L))
	Crossover<-c(rep(0,n))
for(i in 1:n){
if(COS[i]>=Compare[i]){Crossover[i]=1}
if(COS[i]<Compare[i]){Crossover[i]=0}
	}
Total=sum(Crossover)
Pi_Est=(2*n)/Total

cat("For ",n," trials the estimate of pi is: ",Pi_Est)
}
########################################

Buffon_Needle(100)

A SAS implementation of the monte carlo simulation is:

%macro Buffon_Needle(n=);
%let D=2; %let L=2;
 
data BN;
do i=1 to &n;
     angle=rand("Uniform",0,3.14159/2);
     X=rand("Uniform",0,&D/2);
     COS=cos(angle);
     L=&L;
output;
end;
 
data BN; set BN;
if(COS >=(X/(2/L))) then Cross = 1;
else Cross = 0;
 
proc sql;
create table Results as
select (sum(t1.Cross)) as Total_Crosses
from work.BN t1;
quit;
 
data Results; set Results;
n=&n;
pi_est = (2*&n)/Total_Crosses;
 
run;
 
%mend;
 
%Buffon_Needle(n=1000000);

Advertisement

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 )

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