// TP3 - Exercice 1
// Méthode de Monte Carlo - estimation de pi

// On a l'égalité :
// pi = 4 P(U² + V² < 1) = E [4 x 1_{U²+V²<1}]
// où 1_A désigne l'indicatrice de l'évènement A, et U et V sont uniformes sur [0,1].

N = 10000;
U = rand(N,2); // Chaque ligne de U contient un couple de variable
               // uniformes sur [0,1] indépendantes

V = 4 * ((U(:,1).^2 + U(:,2).^2)<1); // Chaque ligne de V contient une
                                     // réalisation de 4 x 1_{U²+V²<1}.
                                     // ( x(:,i) correspond à la ième colone de x )

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


// On affiche l'intervalle de confiance à 95%
disp("Intervalle de confiance à 95% pour pi, 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));


///////////////////////

// Le tracé du graphique

n = [1:N]';
moy_n = cumsum(V)./n; // La i ème case de moy_n contient la moyenne
                      // empirique des i premiers éléments de V.
sigma2_n = cumsum(V.^2)./n - moy_n.^2; // La i ème case de sigma2_n contient
                                   // la variance empirique des i
                                   // premiers éléments de V 

xset('window',0);
clf();
plot2d(n, moy_n, 5); // On trace l'évolution de la moyenne empirique en
                     // fonction de n.

// On trace l'évolution des bornes de l'intervalle de confiance en fonction de n.
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,%pi]);// La valeur théorique
