%% Buffon's needle problem
%% Introduction
% Let us consider the following problem first proposed and solved by the Comte de Buffon in the eighteenth century.
% A needle is randomly thrown on a smooth floor made by parallel strips that have the same width D. Let us call L the length of the needle and suppose that L \leq D. What is the probability that the needle will intersect one of the lines between two strips?
% The problem allows us also to produce a Monte Carlo method to approximate the number \pi.
% Let us consider a generical position of the needle. Let X be the distance from the middle point of the needle to the closest parallel line between strips, and let \theta be the acute angle between the needle and the projected line of length X perpendicular to the strips. If the needle intersects a line, the half needle crossing the line forms a rectangular triangle with the crossed line where the third side is the projected line of length X. This triangle's hypotenuse has length \frac{X}{\cos \theta}. So, for a given position, the needle will intersect one of the parallel lines if and only if \frac{X}{\cos \theta}<\frac{L}{2} \Leftrightarrow X< \frac{L}{2} \cos \theta, otherwise it will completely lie on a strip between two lines.
% X is a random variable that varies uniformly between 0 and \frac{D}{2} while \theta varies uniformly between 0 and \frac{\pi}{2}, and they are independent. Hence the needle crosses a line between two strips with probability
% P\bigg(X<\frac{L}{2} \cos \theta \bigg)= \int_0^{\pi/2} \int_0^{L/2 \cos \theta} \frac{4}{\pi D} dx\, d\theta = \frac{2L}{\pi D} \int_0^{\pi/2} \cos \theta d\theta= \frac{2L}{\pi D}

% Needle drop simulation with L=0.5, D=1
rng('shuffle')
DD=1;
LL=0.5;
if(LL>DD)
    error('ERROR! The parameters must satisfy DD>=LL')
end
Y=rand * DD;
alpha=rand*pi;
figure,hold on
plot([-DD DD],[0 0],'r',[-DD DD],[DD DD],'r')
plot([-cos(alpha)*LL/2 +cos(alpha)*LL/2],[Y-sin(alpha)*LL/2 Y+sin(alpha)*LL/2], 'k','LineWidth',1.4)
ylim([-LL/2 LL/2+DD])
if (alpha<0.5*pi)
    theta=0.5*pi-alpha;
else
    theta=alpha-0.5*pi;
end
if (Y<DD/2)
    X=Y;
else
    X=DD-Y;
end
fprintf(' theta=%g  X=%g  L/2*cos(theta)=%g \n',theta,X,LL/2*cos(theta))
if (Y-sin(alpha)*LL/2<0 || Y+sin(alpha)*LL/2>DD)
    fprintf(' The needle intersect a line')
else
    fprintf(' The needle do not intersect any line')
end

%
%% The Estimate of \pi
% Let us run the Buffon's needle experiment a large number n of times. By the law of large numbers, the ratio of the number of times that the needle crosses a line to the number of trials tends with probability 1 to the probability of a line crossing as n goes to infinity.
% Let us define the random variable N_n that counts the number of line crossings in the first n experiments. On each run we define a Bernoulli variable I_j that takes value 1 if the needle crosses a line in the j-th run and 0 otherwise, for j=1,\ldots,n. Since we are considering independent replications of the same experiment, the random variables I_j are independent and identically distributed and from the previous calculation it follows that P(I_1=1)=\frac{2L}{\pi D}. Thus, the random variable N_n=\sum_{j=1}^n I_j has binomial distribution with parameters n and \frac{2L}{\pi D}.
% The random variable \frac{N_n}{n} has mean \mathbb{E}\bigg(\frac{N_n}{n}\bigg)=\frac{2L}{\pi D} and variance Var\bigg(\frac{N_n}{n}\bigg)=\frac{2L}{n \pi D} \bigg(1- \frac{2L}{\pi D}\bigg) and by the law of large numbers for n large it is highly probable that \frac{N_n}{n}\approx \frac{2L}{\pi D}, hence we estimate \pi as
% \pi \approx \frac{2L\,n}{N_n D}
%
% Let us now produce a Monte Carlo simulation of 1000 runs of the Buffon's needle experiment. We will plot the drawn (\theta,X) as points in the rectangle [0,\pi/2]\times[0,D/2], where we also plot the curve X= \frac{L}{2} \cos \theta. Note that the drawn (\theta,X) that correspond to a line crossing by the needle are those points falling in the area below the curve. The progression of the estimate of \pi during these 1000 replications of the experiment is shown in the successive figure.

% Simulation of 1000 needle throws with L=.5, D=1
niter=1000;
DD=1.;
LL=0.5;
if(LL>DD)
    error('ERROR! The parameters must satisfy DD>=LL')
end
successes=0;
tic
progr1=zeros(niter,1);
figure, hold on
title('The red points represent drawn (\theta,X)')
angle=0:0.01:pi/2;
pl1=plot(angle,cos(angle)*LL/2,'LineWidth', 1.);
ylim([0 DD/2])
xlabel('\theta'), ylabel('X')
for counter=1:niter   
    X=rand * (DD/2);
    theta=rand * pi/2;
    plot(theta,X , 'r.')
    if (X<LL/2*cos(theta))
        successes=successes+1;
    end
    progr1(counter)=2*LL*counter/(successes*DD);
end
hold off
legend(pl1,'X=L cos(\theta) /2')
estim=2*LL*counter/(successes*DD);
fprintf('\n The estimate of pi after %g tosses is %.5f',niter,estim)
figure,
plot(1:niter,progr1,1:niter,pi*ones(niter,1))
legend('estimate evolution','\pi')
xlabel('number of tosses'), ylabel('estimate of pi')
toc

%
% Let us now consider a larger number of iterations of the experiment. The evolution of the estimate of \pi is described in the sequel:

% niter=2x10^8 drops , D=1, L=0.5
niter=2e8;
DD=1.;
LL=.5;
if(LL>DD)
    error('ERROR! The paramenters must satisfy DD>=LL')
end
tic
successes=0;
progr2=zeros(niter/1000,1);
for counter= 1:niter   
    X=rand * (DD/2);
    theta=rand * pi/2;
    if (X<LL/2*cos(theta))
        successes=successes+1;
    end
    if(rem(counter,1000)==0) %saving data each 1000 itarations
        progr2(counter/1000)=2*LL*counter/(successes*DD);
    end
    %if(rem(counter,2e7)==0)
    %    fprintf('\n %g iterations done',counter)
    %end
end
estim=2*LL*counter/(successes*DD);
fprintf('\n The estimate of pi after %g tosses is %.5f',niter,estim)
fprintf('\n The estimate of pi after 10^3 tosses is %.5f',progr2(1))
fprintf('\n The estimate of pi after 10^4 tosses is %.5f',progr2(10))
fprintf('\n The estimate of pi after 10^5 tosses is %.5f',progr2(100))
fprintf('\n The estimate of pi after 10^6 tosses is %.5f',progr2(1000))
fprintf('\n The estimate of pi after 10^7 tosses is %.5f',progr2(10000))
fprintf('\n The estimate of pi after 10^8 tosses is %.5f',progr2(100000))
figure,
plot(1:niter/1000,progr2,1:niter/1000,pi*ones(niter/1000,1))
legend('estimate evolution','\pi')
xlabel('number of tosses x 10^3'), ylabel('estimate')
toc
fprintf('\n The estimate of pi after %g tosses is %.5f',niter,estim)

%
%% Convergence to \pi: efficiency
% Let us now run the experiment an arbitrary number of times for different parameters D and L.
% We want to study the efficiency of this method of estimating \pi. We will make use of the Central Limit Theorem to deduce how many iterations of the needle experiment are needed to produce an estimate of \pi arbitrarly close to the real value of \pi=3.14159265359....
% Let us consider the Bernoulli random variable I_j previously defined. We define a new random variable A_j=\frac{D}{2L} I_j by multiplying I_j by the constant \frac{D}{2L}, for each j=1,2,\ldots,n. Note that A_j takes the value \frac{D}{2L}  if the needle crosses a line in the j-th run, so with probability \frac{2L}{\pi D}, and the value 0 otherwise, with probability 1-\frac{2L}{\pi D}=\frac{\pi D - 2L}{\pi D}. So that the expected value is \mu=\mathbb{E}(A_j)=\frac{1}{\pi} and the variance is \sigma^2=Var(A_j)=\frac{\pi D - 2L}{2L\pi^2}. The A_j are indipendent and identically distributed.
% Defining \overline{S}_n=\frac{1}{n}(A_1+\ldots+A_n) we know that thanks to the Central Limit Theorem, for large n,
% \sqrt n \frac{\overline S_n - \mu }\sigma
% \simeq N(0,1) ,
% where N(0,1) has standard normal distribution. Let us fix a value \varepsilon. We can estimate the probability that \overline S_n differs from \mu less than \varepsilon as
% \mathrm P (|\overline S_n - \mu|\leq \varepsilon)=
% \mathrm P \bigg(\bigg|\frac 1 n (A_1+\ldots+A_n) - \frac 1 \pi \bigg|\leq \varepsilon\bigg)
% \simeq 2 \Phi\Big(\sqrt n \frac \varepsilon \sigma\Big) -1
% where \Phi indicates the cumulative distribution function of N(0,1).
% Note that being \sigma^2 a decreasing function of L\in(0,D], the estimate improves when L grows up to D.
% We can now calculate how many iterations are sufficient to produce \overline S_n = \frac 1 \pi +err, where the error verifies |err|\leq \varepsilon . It holds
% \overline S_n -err =\frac 1 \pi \Leftrightarrow \pi = \frac 1 {\overline S_n - err} = \frac 1 \overline S_n \bigg(1 + \frac {err} {\overline S_n} + o\Big(\frac {err}{\overline S_n} \Big)\bigg)= \frac{1}{\overline{S}_n} + \frac{err}{\overline{S}_n^2}+ o \Big(   \frac{err}{\overline{S}_n^2} \Big)
% so that, if \varepsilon is smaller than 10^{-3 the absolute value of the error in the estimate of \pi is smaller than 10\varepsilon. So, if we want that the estimate \frac 1 {\overline S_n} will differ from \pi of less than 10^{-2}, we must choose \varepsilon=10^{-3}; if we want a difference smaller than 10^{-3} we must choose \varepsilon=10^{-4} and so on.
% Let us evaluate how large must be the number of iterations n to produce an error smaller than an a priori fixed value with high probability. Recalling that \Phi(0.6745)=0.75, \Phi(1.15035)=0.875, \Phi(1.645)=0.95, \Phi(1.96)=0.975 and \Phi(2.58)=0.995, we have that

% number of iterations needed to have error smaller than error_pi
DD=1;
LL=1;
error_pi=1e-3;
if(LL>DD)
    error('ERROR! The parameters must satisfy DD>=LL')
end
sigma_sq=(pi*DD-2*LL)/(2*LL*pi*pi);
eps=error_pi/10;
quant=0.6745;
phi_quant=0.75;
req_iter=floor(sigma_sq*quant^2/eps^2 )  +1;
fprintf('\n D=%g  L=%g ',DD,LL)
fprintf('\nThe absolute value of the error is shorter than %g',error_pi)
fprintf('with probability more than %g after %g iterations',2*phi_quant -1, req_iter)
quant=1.15035;
phi_quant=0.875;
req_iter=floor(sigma_sq*quant^2/eps^2 )  +1;
fprintf('\nThe absolute value of the error is shorter than %g',error_pi)
fprintf('with probability more than %g after %g iterations',2*phi_quant -1, req_iter)
quant=1.645;
phi_quant=0.95;
req_iter=floor(sigma_sq*quant^2/eps^2 )  +1;
fprintf('\nThe absolute value of the error is shorter than %g',error_pi)
fprintf('with probability more than %g after %g iterations',2*phi_quant -1, req_iter)
quant=1.96;
phi_quant=0.975;
req_iter=floor(sigma_sq*quant^2/eps^2 )  +1;
fprintf('\nThe absolute value of the error is shorter than %g',error_pi)
fprintf('with probability more than %g after %g iterations',2*phi_quant -1, req_iter)
quant=2.58;
phi_quant=0.995;
req_iter=floor(sigma_sq*quant^2/eps^2 )  +1;
fprintf('\nThe absolute value of the error is shorter than %g',error_pi)
fprintf('with probability more than %g after %g iterations',2*phi_quant -1, req_iter)

%
% We can now choose the number of iteration and produce the estimate:

niter=5e6;  %10k to 200 mlns
DD=1;
LL=1;
if(LL>DD)
    error('ERROR! The parameters must satisfy DD>=LL')
end
tic
progr3=zeros(niter/1000,1); %save partial data every 1000 iterations in progr
successes=0;
for counter= 1:niter 
    X=rand * (DD/2);
    theta=rand * pi/2;
    if (X<LL/2*cos(theta))
        successes=successes+1;
    end
    if(rem(counter,1000)==0) %save partial data every 1000 iterations
        progr3(counter/1000)=2*LL*counter/(successes*DD);
    end
    if(rem(counter,1e6)==0) %print a message every one million iterations
        fprintf('\n %g iterations done',counter)
    end
end
toc
estim=2*LL*counter/(successes*DD);
figure,
plot(1:niter/1000,progr3,1:niter/1000,pi*ones(niter/1000,1))
xlabel('number of tosses x 10^3'), ylabel('estimate')
legend('estimate evolution','\pi')
fprintf('\n The estimate of pi after %g tosses is %.5f',niter,estim)
fprintf('\nThe absolute value of the error is %g',abs(estim-pi))

%

help rng
help rand

%
%
% Buffon's needle problem - Wikipedia, The Free Encyclopedia
% Monte Carlo method - Wikipedia, the Free Encyclopedia
% Law of large numbers - Wikipedia, The Free Encyclopedia
% Central limit theorem - Wikipedia, The Free Encyclopedia