///////////////////////////////////////////////////////////////////////////////////
//                  ECHANTILLONNAGE D'IMPORTANCE     (SOLUTION EXO2)         //////
///////////////////////////////////////////////////////////////////////////////////


stacksize(100000000);

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

//Nombre de répétitions de l'expérience
nb=1000;

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

//g-Valeur de la sortie g(Y)=g(f(X)), avec f(x)=x et g(x)=(x-2)_+
y=max((x-2),0);

//Les estimateurs empiriques de E(g(Y))=E((Y-2)_+)
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((Y-2)_+) pour Y loi N(0,1), Méthode Monte Carlo']);

end








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


//Calcul de la moyenne E(g(Y))=E((Y-2)_+) par integration 
function y=f1(x), y=(x-2)*exp(-x^2/2)/sqrt(2*%pi),endfunction
I1=intg(2,100,f1)

//Calcul de la variance E((Y-2)_+^2)-E((Y-2)_+)^2 par integration 

function y=f2(x), y=(x-2)^2*exp(-x^2/2)/sqrt(2*%pi),endfunction
I2=intg(2,100,f2)
I=I2-I1^2;


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

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

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




/////////////////////////////////////////////////////////////////////////////////////////////////
///////Histogramme des fluctuations empiriques centrées normées (n=nb=1000) (EXO2,QUESTION 2)////
////////////////////////////////////////////////////////////////////////////////////////////////

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


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

V=(my(n,:)-I1)/sqrt(I);
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")

xtitle(['Histogramme des fluctuations normalisées sur "+string(nb)+" expériences (nombre de simulations par expérience : "+string(n)+") ', 
'Densité gaussienne N(0,1). Méthode Monte Carlo']);

end
end

/////////////////////////////////////////////////////////////////////////////////////////////////
//Echantillonnage d'importance avec une gaussienne centrée en 2, et réduite (EXO2, QUESTION 3)///
/////////////////////////////////////////////////////////////////////////////////////////////////


//Calcul de la moyenne et de la variance E_Q((Y-2)_+^2 (dP/dQ)^2)-E((Y-2)_+)^2 par integration 

function y=f1(x), y=(x-2)*exp(-x^2/2)/sqrt(2*%pi),endfunction
I1=intg(2,100,f1)

function y=f3(x), y=(x-2)^2*exp(-2*x+2)*exp(-x^2/2)/sqrt(2*%pi),endfunction
I3=intg(2,100,f3)
II=I3-I1^2;



xI=grand(n,nb,'nor',2,1);
poids=exp(-xI.^2./2)./(exp(-(xI-2).^2./2));

Iy=max((xI-2),0).*poids;
mIy=cumsum(Iy,'r')./((1:n)'*ones(1,nb));

if nb<20 then

clf(3);xset("window",3);show_window();
plot2d(mIy)
xtitle(['Nombre d''expériences : "+string(nb)+" , Nombre de simulations par expérience : "+string(n)+" ',
'estimation empirique de E((Y-2)_+) pour Y loi N(0,1), Echantillonnage d''importance avec N(2,1)']);

end

/////////////////////////////////////////////////////////////////////////////////////
///////Moyennes et variances empiriques des nb estimateurs "d'importance" simulés////
//////////////////////////////////////////////////////////////////////////////////////
MeanI=mean(Iy,'r');
VarI=variance(Iy,'r');



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

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

plot2d(MeanI);
MeanI1=I1*ones(1,nb);
MeanEmpI=mean(MeanI)*ones(1,nb);
plot2d(MeanI1,style=5);
plot2d(MeanEmpI,style=2);
xtitle(['Courbe des moyennes sur "+string(nb)+" expériences. Nombre de simulations par expérience : "+string(n)+" ', 
'Valeur exacte "+string(I1)+" (rouge), Valeur empirique "+string(mean(MeanI))+" (bleu). Echantillonnage d''importance avec N(2,1)']);

/////////////////////////////////////////////////////////////////////////////////////////////////
///////Histogramme des fluctuations empiriques centrées normées (n=nb=1000) (EXO2,QUESTION 4)////
////                                                            (n=500, nb=100)
////////////////////////////////////////////////////////////////////////////////////////////////

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


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

VI=(mIy(n,:)-I1)/sqrt(II);
VIn=sqrt(n)*VI;
histplot(linspace(min(VIn),max(VIn),20),VIn,style=5);

t = linspace(-5,15,60)'; 
s = exp(-t.^2/2)/sqrt(2*%pi);
  plot2d(t,s,2,"000")
  
  xtitle(['Histogramme des fluctuations normalisées sur "+string(nb)+" expériences (nombre de simulations par expérience : "+string(n)+") ', 
'Densité gaussienne N(0,1). Echantillonnage d''importance avec N(2,1)']);

end
end