///////////////////////////////////////////////////////////////////////////////////////////////////////////
// Simulation de l'aiguille de Buffon (Georges-Louis Leclerc de Buffon (1733)      (SOLUTION EXO1)   //////
///////////////////////////////////////////////////////////////////////////////////////////////////////////

clf();


//Nombre de répétitions de l'expérience de Buffon

I=1000; //Tracé des aiguilles et de l'estimation de pi uniquement pour I=1

//Initialisations pour les vecteurs "erreurs empiriques" et "variables gaussiennes centrées normées"

G=[];
GG=[];

//(Boucle) répétitions des expériences de Buffon

for N=1:I 


//Nombre d'aiguilles simulées
n=10000;

//Simulation des 3 variables uniformes permettant de positionner les abscisses et coordonnées des aiguilles
r=rand(3,n);

//Les coordonnées cartesiennes des aiguilles 
x=[r(1,:);r(1,:)+cos(2*%pi*r(3,:))];
y=[r(2,:);r(2,:)+sin(2*%pi*r(3,:))];
lx=[-1,-1;2,2];
ly=[0,1;0,1];

// rectangle =[xmin,ymin,xmax,ymax]
rectangle=[-1,0,2,1];

///Coloriage des aiguilles en bleu
colors=2*ones(1,n);


//Indices des aiguilles intersectantes + comptage
cuts=find(((y(2,:)>1)|(y(2,:)<0)));
nbcuts=length(cuts);

///Re-coloriage des aiguilles intersecantes en rouge
colors(cuts)=5*ones(cuts);


//Tracé des aiguilles et de l'estimation de 2/pi pour (uniquement) I=1 simulation

if I==1 then

clf(0);xset("window",0);show_window();

//Positionnement de la fenêtre 
plot2d(lx,ly,strf='012',style=[1,1],frameflag=3,axesflag=1,rect=rectangle);

//Tracé des aiguilles
plot2d(x,y,style=colors);

//Titre avec résultat d'estimation
xtitle(["Nombre de lancers:"+string(n)+" Nombre d''intersections:"+string(nbcuts);" Approximation de pi: "+string(2*n/nbcuts)]);


end



///////////////////////////PROGRAMMES NON FOURNIS//////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
////   Estimation par Monte Carlo de 2/pi par les proprotions d'intersections (EXO1, Question 1)///
///////////////////////////////////////////////////////////////////////////////////////////////////

//Selection des aiguilles intersectantes

A=((y(2,:)>1)|(y(2,:)<0));

// 0=True prennent la valeur 1; les autres prennent la valeur 0

A(A==0)=1; 

// Moyennes empiriques au fil du temps
MA=cumsum(A)./(1:n);





///////////////////////////////////////////////////////////////////////////////////////////////////
////   Tracé des estimation empiriques et de la valeur recherchée  (EXO1, Question 2)           ///
///////////////////////////////////////////////////////////////////////////////////////////////////
if I==1 then
clf(1);xset("window",1);show_window();

plot2d(MA);
plot([1:n],2/%pi);
xtitle('Estimation de 2/pi')

end



//Stockage des données à chaque itération [erreurs empiriques centrées et normées]
G=[G,sqrt(n)*(MA(1,n)-2/%pi)/sqrt(((2/%pi)*(1-2/%pi)))];

//Simulation à chaque itération de Gaussiennes iid
GG=[GG,grand(1,1,'nor",0,1)];
end

///////////////////////////////////////////////////////
////   Illustration du TCL  (EXO1, Question 3) ////////
///////////////////////////////////////////////////////

if  I>=100 then

clf(2);xset("window",2);show_window();
histplot(linspace(min(G),max(G),20),G,style=2);
histplot(linspace(min(GG),max(GG),20),GG,style=5);


// Tracé de la courbe gaussienne centrée normée

t = linspace(-5,5,60)'; 
s = exp(-t.^2/2)/sqrt(2*%pi);
  plot2d(t,s,5,"000")
  
xtitle(['Histogrammes des fluctuations normalisées (iid gaussiennes : rouge)' ; 'Fluctuations des "+string(I)+" (simulations = bleu)'])  


//Tracé des fonctions de répartition empiriques
clf(3);xset("window",3);show_window();


//répartition empirique et répartition gaussienne
subplot(121);
Gx=gsort(G,'g','i');
plot2d2(Gx,(1:I)/I,style=3)
y=cdfnor('PQ',Gx,zeros(G),ones(G));
plot2d(Gx,y,style=5);
xtitle('Fonction de répartition "+string(I)+" simulations')  

//répartition empirique et répartition gaussienne
subplot(122);
GGx=gsort(GG,'g','i');
plot2d2(GGx,(1:I)/I,style=3)
y=cdfnor('PQ',GGx,zeros(GG),ones(GG));
plot2d(GGx,y,style=5);
xtitle('Fonction de répartition "+string(I)+" iid gaussiennes')  


//Valeur en u des fonctions de répartition empiriques


u=1.96;
FGu=sum(Gx<=u)/I
FGGu=sum(GGx<=u)/I

//Calcul empirique des proba des intervalles de confiances 


PGd=sum(abs(G)>1.96)/I
PGGd=sum(abs(GG)>1.96)/I





//Inégalité de Mills et estimation de la fonction de répartition gaussienne (sans tables)

clf(4);xset("window",4);show_window();

t= linspace(1,3,60)';
g=exp(-t.^2/2)./sqrt(2*%pi);
tt=t+t^(-1);
  plot2d(t,g./t,style=2)
 plot2d(t,g./tt,style=5)
 plot2d(t,erfc(t./sqrt(2))./2,style=4) 

xtitle('Inégalités de Mills : g(x)/(x+1/x) <P(N(0,1)>x)< g(x)/x') 


/// Encadrement de P(|N(0,1)|< x)

clf(5);xset("window",5);show_window();

plot2d(t,2*(1-g./t)-1,style=2)

plot2d(t,2*(1-g./tt)-1,style=5)

 plot2d(t,2*(1-erfc(t./sqrt(2))./2)-1,style=4) 

xtitle('Inégalités de Mills Encadrement de P(|N(0,1)|< x)') 

end