Exercice 1 : Monte-Carlo

Nous illustrons tout d’abord la méthode de Monte-Carlo sur un exemple très simple.

1. Implémenter une fonction qui génère un échantillon de taille \(10\), \(iid\) selon une loi normale de moyenne \(m=2\) et d’écart-type \(s=3\) et renvoie l’estimation de l’écart-type de cet échantillon.
2. Utiliser la méthode de Monte-Carlo pour estimer l’écart-type de la loi à l’origine de cet échantillon, à partir de multiples réalisations (par exemple 5) de l’estimateur de l’écart-type implémenté ci-dessus. Recommencer avec \(1\,000\) réalisations ? Quel fameux résultat théorique illustre-t-on ici ?

En reprenant l’exemple du cours, programmer un estimation de Monte Carlo du nombre \(\pi\approx 3,1416\).

1. Programmer une fonction roulette_coord() qui prend un seul argument ncases (représentant le nombre de valeurs possible sur la roulette utilisée) valant 35 par défaut, permettant de générer les coordonnées d’un point (entre \(0\) et \(35\)). On pourra utiliser la fonction sample() (consulter l’aide grâce à la commande ?sample). La fonction renverra le vecteur des 2 coordonnées \(x\) et \(y\) ainsi générées.

   roulette_coord <- function(ncases = 35) {
       x <- ...
       y <- ...
       return(c(x, y))
   }

2. En utilisant la formule de la distance entre 2 points, programmer une fonction calculant la distance à l’origine (qui a ici pour coordonnées \((\frac{ncases}{2}, \frac{ncases}{2})\)) : \(d = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}\) et vérifiant si celle-ci est inférieure ou égale au rayon du cercle unité (\(R = \frac{ncases}{2}\)). Cette fonction, qu’on appellera par exemple inside_disk() prendra deux arguments : d’une part le vecteur p des deux coordonnées du point évalué, d’autre part ncases. Elle renverra une valeur booléenne (TRUE ou FALSE) indiquant si le point est effectivement à l’intérieur du disque.

   inside_disk <- function(p, ncases = 35) {
       d <- ...
       return(d <= ...)
   }

3. Comme le ratio entre la surface du disque de rayon \(\frac{ncases}{2}\) et celle du carré de coté \(ncases\) est égal à \(\frac{\pi}{4}\), la probabilité d’échantillonner un point à l’intérieur du disque (sachant qu’il est dans le carré) est de \(\frac{\pi}{4}\). Appuyez-vous sur ce résultat pour programmer une fonction calculant l’estimation de Monte-Carlo de \(\pi\) à partir d’un vecteur booléen de taille \(n\) (le nombre de points échantillonnés) dont chaque élément vaut TRUE si le point correspondant se trouve effectivement à l’intérieur du cercle, et FALSE sinon.

   piMC <- function(in_circle) {
       return(...)
   }

4. À l’aide du code ci-dessous, représenter les données générées puis afficher l’estimation de Monte-Carlo de \(\pi\) correspondante. Faire varier npoints et commenter. (ProTip: essayer ngrid <- 1000 et npoints <- 5000)

   # Taille de la grille (résolution)
   ncases <- 35
   
   # Taille de l'échantillon de Monte Carlo
   npoints <- 200
   
   # Génération des points
   pp <- matrix(NA, ncol = 2, nrow = npoints)
   for (i in 1:nrow(pp)) {
       pp[i, ] <- roulette_coord(ncases)
   }
   
   # Estimateur de Monte-Carlo de Pi
   in_disk <- apply(X = pp, MARGIN = 1, FUN = inside_disk, ncases = ncases)
   piMC(in_disk)
   
   # Grephique on commence par initialiser un plot vide, de la bonne taille avec
   # l'argument `type = 'n'`
   plot(x = pp[, 1], y = pp[, 2], xlim = c(0, ncases), ylim = c(0, ncases), axes = 0,
       xlab = "x", ylab = "y", type = "n")
   ## on gradue les axes x puis y de 0 à ncases
   axis(1, at = c(0:ncases))
   axis(2, at = c(0:ncases))
   ## on ajoute un carré encadrant le dessin
   box()
   ## on trace (en pointillés grâce à l'argument `lty = 3`) la grille sur laquelle
   ## sont échantillonnés nos points
   for (i in 0:ncases) {
       abline(h = i, lty = 3)
       abline(v = i, lty = 3)
   }
   ## on ajoute les points échantillonés
   lines(x = pp[, 1], y = pp[, 2], xlim = c(0, ncases), ylim = c(0, ncases), xlab = "x",
       ylab = "y", type = "p", pch = 16)
   ## on ajoute la représentation du cercle
   x.cercle <- seq(0, ncases, by = 0.1)
   y.cercle <- sqrt((ncases/2)^2 - (x.cercle - ncases/2)^2)
   lines(x.cercle, y = y.cercle + ncases/2, col = "red")
   lines(x.cercle, y = -y.cercle + ncases/2, col = "red")
   ## enfin on colorie en rouge les points échantillonnés à l'intérieur du cercle
   lines(x = pp[in_disk, 1], y = pp[in_disk, 2], xlim = c(0, ncases), ylim = c(0, ncases),
       xlab = "x", ylab = "y", type = "p", pch = 16, col = "red", cex = 0.7)