// TP3 - Exercice 2
// Méthode de Monte Carlo - estimation de \int_0^a e^{-x^2}dx
// La valeur théorique est de 0.845903... pour a=sqrt(2), et de 0.882081... pour a=2.

N = 10000;
a = sqrt(2);
// a = 2;


///////////// Méthode 1 : avec une variable de loi uniforme /////////////

// On a l'égalité :
// \int_0^a e^{-x^2}dx = a \int_0^1 e^{-(ax)^2}dxa = a E[exp(-(aU)^2)]
// où U est uniforme sur [0,1].


U = rand(N,1);
X = a * exp(-(a*U).^2);

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

// L'intervalle de confiance
disp("Méthode 1 :");
disp("Intervalle de confiance à 95% pour I , 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 "+string(sqrt(%pi)/2*erf(a)));
disp("");

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

// On trace
xset('window',1);
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], sqrt(%pi)/2*erf(a)*[1 1]); // La valeur théorique



///////////// Méthode 2 : avec une variable Gaussienne /////////////

// On a l'égalité :
// \int_0^a e^{-x^2}dx = sqrt(pi) * \int_0^(sqrt(2)*a) exp(-x^2)dx/sqrt(2pi) = sqrt(pi)P(0<G<sqrt(2)a) = sqrt(pi)*E[1_[0,sqrt(2)a](G)]
// où G est de loi normale N(0,1)


G = grand(N,1,'nor',0,1);
X = sqrt(%pi)*(G>0).*(G<sqrt(2)*a);

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

// L'intervalle de confiance
disp("Méthode 2 :");
disp("Intervalle de confiance à 95% pour I , 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 "+string(sqrt(%pi)/2*erf(a)));
disp("");

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

// On trace
xset('window',2);
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], sqrt(%pi)/2*erf(a)*[1 1]); // La valeur théorique



///////////// Méthode 3 : autre manière avec une variable Gaussienne /////////////

// On a l'égalité :
// \int_0^a e^{-x^2}dx = 1/2*\int_-a^a e^{-x^2}dx =  sqrt(pi)/2*E[1_[-sqrt(2)a,sqrt(2)a](G)]
// où G est de loi normale N(0,1)


G = grand(N,1,'nor',0,1);
X = sqrt(%pi)/2*(abs(G)<sqrt(2)*a);

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

// L'intervalle de confiance
disp("Méthode 3 :");
disp("Intervalle de confiance à 95% pour I , 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 "+string(sqrt(%pi)/2*erf(a)));
disp("");

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

// On trace
xset('window',3);
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], sqrt(%pi)/2*erf(a)*[1 1]); // La valeur théorique



disp("La meilleure méthode est celle pour laquelle la largeur de l intervalle est la plus petite. Ici, il s agit de la méthode 3. La méthode 2 est la plus mauvaise.");
