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

Advertisements