////////////////////////////////////////////////////////////////////////////////////////////
//   STRATIFICATION - ESTIMATION de E(exp(X)) pour X uniforme [-1,1]  (SOLUTION EXO3)  //////
///////////////////////////////////////////////////////////////////////////////////////////

stacksize(268435454);

//Nombre de simulation (à chaque expérience) CHOISIR UN NOMBRE n = 4 m 

n=1000;

//Nombre de répétitions de l'expérience = nb
//Tracé des estimations empiriques sur nb expériences uniquement pour un nombre de simulation 
//inférieur à 50

nb=1000;

//Simulation des variables uniformes sur [-1,1]
x=grand(n,nb,'unf',-1,1);

//g-Valeur de la sortie g(Y)=g(f(X)), avec f(x)=exp(x) et g(x)=x
y=exp(x);

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

//Tracé des estimations empiriques sur nb expériences uniquement pour un nombre de simulation 
// inférieur à 50

if nb <= 50 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 [-1,1]. Méthode de Monte Carlo']);

end

///////////////////////////PROGRAMMES NON FOURNIS/////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////
///////Moyennes et variances empiriques des nb estimateurs simulés (EXO3,QUESTION 1)////
////////////////////////////////////////////////////////////////////////////////////////
Mean=mean(y,'r');
Var=variance(y,'r');
MeanVar=mean(Var);


///TRACE des courbes des moyennes////////

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

plot2d(Mean);
MeanE=mean(Mean)*ones(1,nb);
MeanI=sinh(1)*ones(1,nb);
plot2d(MeanI,style=5);
plot2d(MeanE,style=2);
xtitle(['Courbe des moyennes sur "+string(nb)+" expériences. Nombre de simulations par expérience : "+string(n)+" ', 
'Valeur exacte "+string(sinh(1))+" (rouge), Valeur empirique "+string(mean(Mean))+" (bleu). Méthode Monte Carlo']);



//////////////////////////////////////////////////////////////////////////////////////////////
////     VARIABLES ANTITHETIQUES  f(u)=e^u et f(u)=e^(-u)  sur [-1,1]   (EXO3,QUESTION 2)  ///
//////////////////////////////////////////////////////////////////////////////////////////////

m=n/2;
x_anti=grand(m,nb,'unf',-1,1);


// Couplage des simulation selon la méthode des variables antithétques
y_anti1=exp(x_anti);
y_anti2=exp(-x_anti);

//Les estimateurs empiriques de E(g(Y))=E(e^Y); et leur moyenne sur nb expériences

y_anti=(y_anti1+y_anti2)/2;
my_anti=cumsum(y_anti,'r')./((1:m)'*ones(1,nb));
Mmy_anti=mean(my_anti,'r');

//Variances des différentes expériences, et variance moyennes
Var_anti=variance(y_anti,'r');
MVar_anti=mean(Var_anti);


///TRACE des courbes des moyennes

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

plot2d(Mmy_anti);
MeanI=sinh(1)*ones(1,nb);
Mean_Anti=mean(Mmy_anti)*ones(1,nb);
plot2d(Mean_Anti,style=2);
plot2d(MeanI,style=5);
xtitle(['Courbe des moyennes sur "+string(nb)+" expériences. Nombre de simulations par expérience : "+string(n)+" ', 
'Valeur exacte "+string(sinh(1))+" (rouge), Valeur empirique "+string(mean(Mmy_anti))+" (bleu) Variables antithétques']);




////////////////////////////////////////////////////////////////////
//Simulation stratification proportionnelle (EXO3,QUESTION 3)//////
///////////////////////////////////////////////////////////////////
n1=n/2;
n2=n1;


x1=grand(n1,nb,'unf',-1,0);
x2=grand(n2,nb,'unf',0,1);

y1=exp(x1);
y2=exp(x2);
z=(y1+y2);


//Les estimateurs empiriques de E(g(Y))=E(e^Y)

mz=cumsum(z,'r')./((1:n1)'*ones(1,nb));


//Tracé des estimations empiriques sur nb expériences uniquement pour un nombre de simulation 
// inférieur à 50

if nb <= 50 then

clf(3);xset("window",3);show_window();
plot2d(mz)
xtitle(['Nombre d''expériences : "+string(nb)+" , Nombre de simulations par expérience : "+string(n)+" ',
'Estimation empirique de E(exp(X)) pour X uniforme [-1,1], Stratification proportionnelle sur [-1,0]U[0,1]']);

end


Mean1=mean(z,'r')/2;
Var1=variance(z,'r')/2;
MeanVar1=mean(Var1);


///TRACE des courbes des moyennes////////

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

plot2d(Mean1);
MeanI=sinh(1)*ones(1,nb);
MeanE1=mean(Mean1)*ones(1,nb);
plot2d(MeanE1,style=2);
plot2d(MeanI,style=5);
xtitle(['Courbe des moyennes sur "+string(nb)+" expériences (nombre de simulations par expérience) : "+string(n)+" ', 
'Valeur exacte "+string(sinh(1))+" (rouge), Valeur empirique "+string(mean(Mean))+" (bleu). Stratification proportionnelle sur [-1,0]U[0,1]']);



//////////////////////////////////////////////////////////////////////
//Simulation stratification non proportionnelle (EXO3,QUESTION 4)//////
///////////////////////////////////////////////////////////////////////

rn=1/20;

n1=n*rn;
n2=n*(1-rn);



xx1=grand(n1,nb,'unf',-1,0);
xx2=grand(n2,nb,'unf',0,1);

yy1=exp(xx1);
yy2=exp(xx2);




//Les estimateurs empiriques de E(g(Y))=E(e^Y)

myy1=cumsum(yy1,'r')./((1:n1)'*ones(1,nb));
myy2=cumsum(yy2,'r')./((1:n2)'*ones(1,nb));


//Tracé des estimations empiriques sur nb expériences uniquement pour un nombre de simulation 
// inférieur à 50

if nb <= 50 then

clf(5);xset("window",5);show_window();
plot2d(myy1)
plot2d(myy2)
xtitle(['Nombre d''expériences : "+string(nb)+" , Nombre de simulations par expérience : "+string(n)+" ',
'Estimation empirique de E(exp(X)) pour X uniforme [-1,1], Stratification "+string(n1)+" sur [-1,0] et "+string(n2)+" sur U[0,1]']);

end


Mean2=(mean(yy1,'r')+mean(yy2,'r'))/2;
Var2=(variance(yy1,'r')/rn+variance(yy2,'r')/(1-rn))/4;
MeanVar2=mean(Var2);


///TRACE des courbes des moyennes////////

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

plot2d(Mean2);
MeanI=sinh(1)*ones(1,nb);
MeanE2=mean(Mean2)*ones(1,nb);
plot2d(MeanE2,style=2);
plot2d(MeanI,style=5);

xtitle(['Courbe des moyennes sur "+string(nb)+" expériences. Nombre de simulations par expérience : "+string(n)+" ', 
'Valeur exacte "+string(sinh(1))+" (rouge), Valeur empirique "+string(mean(Mean2))+" (bleu). Stratification "+string(n1)+" sur [-1,0] et "+string(n2)+" sur U[0,1]']);


///////////////////////////////////////////////////////////////////////
//Simulation stratification optimale (EXO3,QUESTION 5)           //////
///////////////////////////////////////////////////////////////////////


sigma1=exp(-1)*(2-3/(2*exp(1)))-1/2;
sigma2=2*exp(1)*(1-exp(1)/4)-3/2;
r=round(10*(sqrt(sigma1)/(sqrt(sigma1)+sqrt(sigma2))))/10;

n1=n*r;
n2=n*(1-r);



xo1=grand(n1,nb,'unf',-1,0);
xo2=grand(n2,nb,'unf',0,1);

yo1=exp(xo1);
yo2=exp(xo2);



//Les estimateurs empiriques de E(g(Y))=E(e^Y)

myo1=cumsum(yo1,'r')./((1:n1)'*ones(1,nb));
myo2=cumsum(yo2,'r')./((1:n2)'*ones(1,nb));


//Tracé des estimations empiriques sur nb expériences uniquement pour un nombre de simulation 
// inférieur à 50

if nb <= 50 then

clf(7);xset("window",7);show_window();
plot2d(myo1)
plot2d(myo2)

end


Mean3=(mean(yo1,'r')+mean(yo2,'r'))/2;
Var3=(variance(yo1,'r')/r+variance(yo2,'r')/(1-r))/4;
MeanVar3=mean(Var3);
xtitle(['Nombre d''expériences : "+string(nb)+" , Nombre de simulations par expérience : "+string(n)+" ',
'Estimation empirique de E(exp(X)) pour X uniforme [-1,1], Stratification optimale "+string(n1)+" sur [-1,0] et "+string(n2)+" sur U[0,1]']);



///TRACE des courbes des moyennes////////

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

plot2d(Mean3);
MeanI=sinh(1)*ones(1,nb);
MeanE3=mean(Mean3)*ones(1,nb);
plot2d(MeanE3,style=2);
plot2d(MeanI,style=5);
xtitle(['Courbe des moyennes sur "+string(nb)+" expériences. Nombre de simulations par expérience : "+string(n)+" ', 
'Valeur exacte "+string(sinh(1))+" (rouge), Valeur empirique "+string(mean(Mean3))+" (bleu). Stratification optimale "+string(n1)+" sur [-1,0] et "+string(n2)+" sur U[0,1]']);




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


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


MeanVart=mean(MeanVar)*ones(1,nb);  ///MONTE CARLO
MeanVart1=mean(MeanVar1)*ones(1,nb); //STRATIF. PROPORTIONNELLE
MeanVart2=mean(MeanVar2)*ones(1,nb); //STRATIF. NON PROPORTIONNELLE
MeanVart3=mean(MeanVar3)*ones(1,nb); //STRATIF. OPTIMALE
MeanVar_anti=mean(MVar_anti)*ones(1,nb); //VARIABLES ANTITHETIQUE


plot2d(MeanVart,style=2);
plot2d(Var,style=1);

plot2d(MeanVart1,style=3);
plot2d(Var1,style=5);

plot2d(MeanVart2,style=4);
plot2d(Var2,style=7);

plot2d(MeanVart3,style=6);
plot2d(Var3,style=9);

plot2d(MeanVar_anti,style=5);
plot2d(Var_anti,style=1);


xtitle(['Comparaison des trois variances ("+string(nb)+" expériences de "+string(n)+" simulations, Monte Carlo [noir])'; 
'(Statif proportionnelle [rouge] non proportionnelle[jaune], stratif optimale [blue], antithétique [noir, centre rouge])']);