// TP3 - Exercice 3
// Méthode de Monte Carlo - estimation de la moyenne de E[e^(X^2)],
// où X suit la loi de densité f(x)=2/3(x*1_[0,1] + 1_[1,2])
// La valeur théorique est 10.56608...


// On simule X par la méthode du rejet :
// f(x) < 4/3g(x),  où  g=1_[0,2]/2 est la densité de la loi uniforme sur [0,1]


N = 10000; // Nombre de simulations

X = []
while size(X,'c')<N
  x = 2*rand(); // x est uniforme sur [0,2]
  if x>1 | rand()<x // Dans la méthode du rejet, l'inégalité 4/3*g(x)*U<f(x) se récrit x>1 ou U<x 
	X = [X, x];
  end
end

T = exp(X.^2)

moy = sum(T)/N; // La moyenne empirique
sigma2 = sum(T.^2)/N - moy^2; // La variance empirique

// L'intervalle de confiance
disp("Intervalle de confiance à 95% pour E[e^(X^2)], avec "+string(N)+" simulations :");
disp("[ " + string(moy-1.96*sqrt(sigma2/N)) + " ; " + string(moy+1.96*sqrt(sigma2/N)) + " ]");
disp("L intervalle est centré en "+string(moy)+" et a une largeur de "+string(2*1.96*sqrt(sigma2/N)));
disp("La valeur exacte est de 10.566078.");

n = [1:N];
moy_n = cumsum(T)./n; // Moyenne empirique en fonction de n
sigma2_n = cumsum(T.^2)./n - moy_n.^2; // Variance empirique en fonction de n  

// On trace l'intervalle de confiance en fonction du nombre de simulations
xset('window',0);
clf();
plot2d(n, moy_n, 5);
plot2d(n, moy_n-1.96*sqrt(sigma2_n./n), 13);
plot2d(n, moy_n+1.96*sqrt(sigma2_n./n), 13);
plot2d([0,N],10.566078*[1 1])
