// TP3 - Exercice 4
// Méthode de Monte Carlo - estimation de
// I = \int_1^\infty dx/(1+x^3) = \int_1^\infty (x^3/2(1+x^3)) 2dx/x^3,
// Autrement dit I = E[f(X)] où X suit la loi de Pareto de paramètre 2 (de densité 2x^-3) et f(x) = x^3*(2+2x^3)^-1
// La valeur exacte est (pi/sqrt(3)-log(2))/3 = 0.37355...
//
// Par le changement de variable y=1/x, on a aussi
// I = \int_0^1 ydy/(1+y^3), de sorte que
// I = E[g(Y)] avev Y uniforme sur [0,1] et g(x)=x/(1+x^3).



N = 10000; // Nombre de simulations



//////////// Méthode 1 : loi de Pareto ////////////

X = rand(N,1).^(-1/2) // On simule une loi de Pareto de paramètre 2.
T = X.^3 ./ (2+2*X.^3)


moy = sum(T)/N; // La moyenne empirique
sigma2 = sum(T.^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((%pi/sqrt(3)-log(2))/3));
disp("");

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',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],(%pi/sqrt(3)-log(2))/3*[1 1])



//////////// Méthode 2 : loi uniforme ////////////

X = rand(N,1) // On simule une loi de Pareto de paramètre 2.
T = X ./ (1+X.^3)


moy = sum(T)/N; // La moyenne empirique
sigma2 = sum(T.^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((%pi/sqrt(3)-log(2))/3));
disp("");

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',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],(%pi/sqrt(3)-log(2))/3*[1 1])



disp("La meilleure méthode est celle qui mène à l intervalle le moins large. Ici, il s agit de la premiere.")
