////////////////////////////////////////////////////////////////////////////////////////////////////
//   VARIABLES DE CONTRÔLE - ESTIMATION de E(exp(X)) pour X uniforme [0,1]   (SOLUTION EXO4)  //////
////////////////////////////////////////////////////////////////////////////////////////////////////

stacksize(100000000);

//Nombre de simulation (à chaque expérience)
n=500;

//Nombre de répétitions de l'expérience. 
//Tracé des fluctuations uniquement pour n>500, et nb>1000
//Tracé des différentes expériences uniquement pour nb <50
nb=50;

//Simulation des gaussiennes centrées réduites
x=grand(n,nb,'unf',0,1);

//g-Valeur de la sortie g(Y)=g(f(X)), avec f(x)=e^x, f_r(x)=1+x et g(x)=x
y=exp(x);
yc=(3/2)+exp(x)-(1+x);


//Les estimateurs empiriques de E(g(Y))=E(e^X)
my=cumsum(y,'r')./((1:n)'*ones(1,nb));
myc=cumsum(yc,'r')./((1:n)'*ones(1,nb));

if nb<51 then

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

xtitle(['Nombre d''expériences : "+string(nb)+" , Nombre de simulations par expérience : "+string(n)+" ',
'estimation empirique de E(exp(X)) pour X uniforme [0,1]. Méthode de Monte Carlo']);


clf(1);xset("window",1);show_window();
plot2d(myc)

xtitle(['Nombre d''expériences : "+string(nb)+" , Nombre de simulations par expérience : "+string(n)+" ',
'estimation empirique de E(exp(X)) pour X uniforme [0,1]. Variable de controle E(1+X)=3/2']);


end


Mean=mean(y,'r');
Meanc=mean(yc,'r');
Var=variance(y,'r');
Varc=variance(yc,'r');


/////////////////////////////////////////////////////////////////////////////////////////////////
///             TRACE des courbes des moyennes empiriques (EXO4,QUESTION 2)                 ////
////////////////////////////////////////////////////////////////////////////////////////////////

clf(2);xset("window",2);show_window();
xtitle(['(Monte Carlo, Var. contrôle,Exact)=(rouge,blue,vert)',
'Moyennes de"+string(nb)+" expériences ("+string(n)+" simul./exp)']);

subplot(121);
plot2d(Mean,style=5);
MeanEmp=mean(Mean)*ones(1,nb);
MeanExact=(exp(1)-1)*ones(1,nb);
plot2d(MeanEmp,style=5);
plot2d(MeanExact,style=3);



subplot(122);
plot2d(Meanc,style=2);
MeanEmpc=mean(Meanc)*ones(1,nb);
MeanExact=(exp(1)-1)*ones(1,nb);
plot2d(MeanEmpc,style=2);
plot2d(MeanExact,style=3);



///////////////////////////PROGRAMMES NON FOURNIS/////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////
///////Histogramme des fluctuations empiriques centrées normées (n=nb=1000) (EXO4,QUESTION 2)/////
//////                                                          (n,nb)=(500,1000)           //////
//////////////////////////////////////////////////////////////////////////////////////////////////

moy=exp(1)-1;
sigma=2*exp(1)-exp(2)/2-3/2;
sigmac=3*exp(1)-exp(2)/2-53/12;

if n>=500 then
    if nb>=1000 then


clf(3);xset("window",3);show_window();
xtitle(['Monte Carlo (gauche)  Var. contrôle (droire)',
'Histo. des fluctuations norm. ("+string(n)+" sim./exp "+string(nb)+")']);

subplot(121);
V=(my(n,:)-moy)/sqrt(sigma);
Vn=sqrt(n)*V;
histplot(linspace(min(Vn),max(Vn),20),Vn,style=5);

t = linspace(-5,15,60)'; 
s = exp(-t.^2/2)/sqrt(2*%pi);
  plot2d(t,s,2,"000")
  
  
subplot(122);
Vc=(myc(n,:)-moy)/sqrt(sigmac);
Vcn=sqrt(n)*Vc;
histplot(linspace(min(Vcn),max(Vcn),20),Vcn,style=5);

t = linspace(-5,15,60)'; 
s = exp(-t.^2/2)/sqrt(2*%pi);
  plot2d(t,s,2,"000")  

end
end

////////////////////////////////////////////////////////////////////
////COMPARAISON DES VARIANCES       (même nombre de simulation)/////
////////////////////////////////////////////////////////////////////


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


xtitle(['(Monte Carlo, Var. contrôle,Exact)=(rouge,blue,vert)',
'Variances de"+string(nb)+" expériences ("+string(n)+" simul./exp)']);

MeanVar=mean(Var)*ones(1,nb);  ///MONTE CARLO
MeanVarc=mean(Varc)*ones(1,nb); //VARIABLE DE CONTRÔLE

subplot(121);
plot2d(Var,style=5);
plot2d(MeanVar,style=3);
subplot(122);
plot2d(Varc,style=2);
plot2d(MeanVarc,style=3);