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!\

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);

### Like this:

Like Loading...

*Related*