3 Simulations de Monte-Carlo

Il est souvent utilse de faire des simulations de Monte-Carlo afin d’étudier les propriétés des estimateurs que nous rencontrerons, et de “visualiser” leurs comprtements. Ces propriétés (biais, convergence, efficacité etc.) sont des propriétés statistques qui ont trait à la distribution des estimateurs, vus comme une variables aléatoire. Une simulation de Monte-Carlo consiste en quelque sortes à tirer plein de valeurs de ces variables aléatoires afin d’étudier leurs distribution.

3.1 Un premier exemple

Votre cours de statistrique et de probabilité vous a (normalement) appris que la moyenne d’un échantillon de taille \(n\) issu d’une loi Normale \(N(\mu,\sigma^2)\) de moyenne \(\mu\) et d’écart-type \(\sigma\) est une variable aléatoire suivant une loi normale de moyenne \(\mu\) et d’écart-type \(\frac{\sigma}{\sqrt{n}}\). Utilisons R pour le vérifier.

Commençons par calculer la moyenne d’un échantillon de 100 tirages d’une \(N(0,1)\)

n <- 100
mu <- 0
sigma <- 1
x <- rnorm(n,mean=mu,sd=sigma)
mean(x)

## [1] -0.03247968

On voit que la moyenne des 100 tirages n’est pas strictement égale à 0, et on ne voit pas bien comment juger de l’écart-type ou de la loi de probabilité de cette moyenne. Si on lance le code une seconde fois, le résultat va d’ailleurs changer :

n <- 100
mu <- 0
sigma <- 1
x <- rnorm(n,mean=mu,sd=sigma)
mean(x)

## [1] -0.1270524

La raison est que la théorie nous donne les caractéristiques de la distribution de cette moyenne au travers d’un grand nombre d’échantillons de 1000 tirages d’une \(N(0,1)\). Nous allons donc répéter le code ci-dessous un grand nombre de fois (\(K\) fois), noter à chaque fois la moyenne obtenue, et étudier la distribution de ces \(K\) moyennes.

Une première façon de faire consiste à utiliser une boucle for :

K <- 10000 # On va faire 10000 réplications, on aura donc 10000 moyennes
n <- 100
mu <- 0
sigma <- 1
moyennes <- c() # On crée un vecteur vide qui contiendra les K moyennes
for (i in 1:K) { # on initialise la boucle
  x <- rnorm(n,mean=mu,sd=sigma) # on tire l'ééchantillon
  moyennes[i] <- mean(x) # on stocke la moyenne issue de la ième réplication à
                          # la ième position du vecteur "moyennes"
}
moyennes[1:10] # on affiche les 10 premières moyennes

##  [1]  0.04282867  0.21528147 -0.09064180 -0.02316018  0.01417983  0.01698663
##  [7]  0.17348736 -0.01235884  0.23749812 -0.11105861

R propose toutefois une fonction dédiée à ce genre de répétition : replicate(). Elle prend deux arguments principaux : le nombre de réplications, et l’expression à répéter. Elle renvoie directement le vecteur des résultats, ce qui nous évite d’écrire la boucle et de gérer nous-mêmes le vecteur moyennes.

K <- 10000
n <- 100
mu <- 0
sigma <- 1
# À chaque réplication, on tire un échantillon et on en calcule la moyenne
moyennes <- replicate(K, mean(rnorm(n, mean = mu, sd = sigma)))
moyennes[1:10] # on affiche les 10 premières moyennes

##  [1] -0.08714515  0.08825767 -0.26370478 -0.08185618 -0.09904880 -0.06632692
##  [7] -0.15473979  0.17106159 -0.09142349 -0.10830222

mean(moyennes) # moyenne des moyennes (0 en théorie)

## [1] 7.841688e-05

sd(moyennes) # écart-type des moyennes (1/sqrt(n) en théorie)

## [1] 0.09971326

1/sqrt(n) # valeur de 1/sqrt(n)

## [1] 0.1

hist(moyennes,prob=TRUE) # on trace l'histogramme des moyennes
lines(density(moyennes),col="red") # on y ajoute le tracé de la densité des moyennes

Ceux qui ont suivi le tutoriel “Visualize Data” peuvent utiliser les outils graphiques du package ggplot2

library(ggplot2)
ggplot(mapping=aes(x=moyennes)) +
         geom_histogram(aes(y=..density..),fill="grey",color="black") +
          geom_density(fill="blue",alpha=0.2)

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

NB : le y=..density.. indique à geom_histogram d’utiliser la densité en ordonnées, au lieu de la fréquence. Ça permet d’avoir l’histogramme et la courbe des densité à la même échelle.

On constate que 1. La moyenne des K=10000 moyennes est très proche de la moyenne théorique (0) 2. L’écart-type des K moyennes est très proche de l’écart-type théorique (0.1) 3. La distribution est proche de celle d’une loi normale

3.2 Généralisation : faire varier les paramètres

L’exemple précédent a montré que l’on pouvait assez facilement effectuer une simulation basique avec du code assez simple. Néanmoins, ce code ne permet de simuler qu’un seul scénario (valeurs de \(n\), de \(\mu\) et de \(\sigma\)) à la fois. Si on souhaite voir ce qu’il se passe lorsqu’on fait varier les paramètres de la simulation, il faudrait copier-coller le code de nombreuses fois et modifer à chaque fois les paramètres. De plus, il faudrait faire attention à sauvegarder les résultats dans des vecteurs différents à chaque fois.

Afin de s’éviter ces désagréments, nous allons automatiser ces tâches en combinant deux outils très répandus : replicate(), que nous venons de voir, et les fonctions map() / pmap() du package purrr (inclus dans le tidyverse), qui servent à appliquer une même fonction à toute une série de valeurs.

On commence par charger le tidyverse (voir les tutoriels)

library(tidyverse) # charge plein d'outils utiles, dont ggplot2, dplyr, tidyr, purrr...

On définit ensuite une fonction qui va effectuer le tirage de l’échantillon aléatoire, calculer la moyenne, et retourner le résultat :

simul_moyenne <- function(n,mu,sigma) { # Notre fonction s'appelle "simul_moyenne" et prend les
                                        # arguments n, mu et sigma
  tirages <- rnorm(n,mean=mu,sd=sigma) # on effectue le tirage aléatoire avec les valeurs données en argument
  mean(tirages) # on retourne la moyenne
}

On définit ensuite l’ensemble des “scénarios” de simulations, c’est à dire toutes les combinaisons des valeurs des paramètres que l’on veut faire varier. Ici on veut deux valeurs pour la taille d’échantillon, deux valeurs pour l’espérance, et trois pour l’écart-type. La fonction expand_grid() (du package tidyr, inclus dans le tidyverse) construit pour nous la data frame de toutes les combinaisons possibles :

scenarios <- expand_grid(n = c(10, 100),
                         mu = c(0, 5),
                         sigma = c(1, 2, 4))
scenarios

## # A tibble: 12 × 3
##        n    mu sigma
##    <dbl> <dbl> <dbl>
##  1    10     0     1
##  2    10     0     2
##  3    10     0     4
##  4    10     5     1
##  5    10     5     2
##  6    10     5     4
##  7   100     0     1
##  8   100     0     2
##  9   100     0     4
## 10   100     5     1
## 11   100     5     2
## 12   100     5     4

On obtient bien \(2 \times 2 \times 3 = 12\) scénarios, chacun correspondant à une ligne.

Il ne reste plus qu’à lancer, pour chaque scénario, replicate() afin d’obtenir les 10000 moyennes correspondantes. C’est le rôle de pmap() : cette fonction parcourt les lignes de notre data frame de scénarios et, pour chacune, appelle une fonction en lui passant les valeurs de n, mu et sigma de la ligne. Comme chaque scénario produit un vecteur de 10000 moyennes, on stocke ce vecteur dans une nouvelle colonne (une colonne-liste), que l’on “déplie” ensuite avec unnest() pour retrouver une data frame classique avec une ligne par réplication.

K <- 10000 # nombre de réplications par scénario

data_resultats <- scenarios %>%
  mutate(moyenne = pmap(list(n, mu, sigma), # pour chaque ligne (n, mu, sigma)...
                        \(n, mu, sigma) replicate(K, simul_moyenne(n, mu, sigma)))) %>% # ...on fait K réplications
  unnest(moyenne) # on déplie : une ligne par réplication

NB : la notation \(n, mu, sigma) ... est une façon abrégée d’écrire function(n, mu, sigma) ... (disponible depuis R 4.1). Les deux écritures sont strictement équivalentes ; on aurait donc aussi pu écrire function(n, mu, sigma) replicate(K, simul_moyenne(n, mu, sigma)).

Regardons ce que l’objet “data_resultats” contient. la fonction head() permet de lister les premières lignes d’une data frame, afin de vérifier rapidement son contenu.

head(data_resultats)

## # A tibble: 6 × 4
##       n    mu sigma moyenne
##   <dbl> <dbl> <dbl>   <dbl>
## 1    10     0     1 -0.145 
## 2    10     0     1 -0.0995
## 3    10     0     1  0.204 
## 4    10     0     1 -0.219 
## 5    10     0     1  0.0430
## 6    10     0     1 -0.606

dim(data_resultats) # 12 scénarios x 10000 réplications = 120000 lignes

## [1] 120000      4

On obtient une data frame bien pratique : une colonne par paramètre (n, mu, sigma) et une colonne moyenne contenant le résultat de chaque réplication.

On va maintenant calculer les moyennes et écart-types des 10000 moyennes calculées pour chaque scénario de simulation. On fait appel aux outils de “pipe” (%>%) et de groupe (group_by()) expliqués dans le tutoriel “Work With Data”. On en profite pour ajouter une variable “ecty_theor” donnant l’écart-type théorique de la distribution des moyennes

data_resultats %>%
mutate(ecty_theor=sigma/sqrt(n)) %>%
group_by(n,mu,sigma) %>%
  summarise(moy=mean(moyenne),ecty=sd(moyenne),ecty_theor=mean(ecty_theor))

## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by n, mu, and sigma.
## ℹ Output is grouped by n and mu.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(n, mu, sigma))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

## # A tibble: 12 × 6
## # Groups:   n, mu [4]
##        n    mu sigma        moy   ecty ecty_theor
##    <dbl> <dbl> <dbl>      <dbl>  <dbl>      <dbl>
##  1    10     0     1 -0.000172  0.314       0.316
##  2    10     0     2  0.00683   0.625       0.632
##  3    10     0     4  0.00675   1.26        1.26 
##  4    10     5     1  5.00      0.314       0.316
##  5    10     5     2  5.00      0.628       0.632
##  6    10     5     4  4.98      1.25        1.26 
##  7   100     0     1 -0.000298  0.0996      0.1  
##  8   100     0     2 -0.0000693 0.198       0.2  
##  9   100     0     4  0.00808   0.400       0.4  
## 10   100     5     1  5.00      0.0995      0.1  
## 11   100     5     2  5.00      0.201       0.2  
## 12   100     5     4  5.00      0.399       0.4

On constate que la moyenne des moyenne (colonne “moy”) est très proche de la moyenne théorique (“mu”), de même que l’écart-type de la distribution des moyennes (“ecty”) est très proche de l’écart-type théorique. Nos observations semblent bien coller avec la théorie.

On complète l’exercice en faisant un graphique de la fonction de densité estimées de nos résultats, pour chaque valeur de \(n\), \(\mu\) et \(\sigma\).

ggplot(data=data_resultats,mapping=aes(x=moyenne,group=sigma,color=factor(sigma))) +
  geom_density() +
  facet_wrap(mu ~ n,labeller=label_both)

3.3 Avec une régression par MCO

R permet de faire une régression linéaire par MCO avec la commande lm()(voir le chapitre sur la répression par MCO) ci-dessous

On va charger une base de données préinstallée avec R : les données “mtcars”

data(mtcars)
head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Effectuons une régression linaire de, par exemple, “mpg” sur “cyl” et “disp”

lm(data=mtcars,mpg ~ cyl + disp)

## 
## Call:
## lm(formula = mpg ~ cyl + disp, data = mtcars)
## 
## Coefficients:
## (Intercept)          cyl         disp  
##    34.66099     -1.58728     -0.02058

La sortie est assez minimale. Stockons cette régression dans un objet que nous appellerons “ma_regression” et faisons un summary() de ce dernier

ma_regression <- lm(data=mtcars,mpg ~ cyl + disp)
summary(ma_regression)

## 
## Call:
## lm(formula = mpg ~ cyl + disp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4213 -2.1722 -0.6362  1.1899  7.0516 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.66099    2.54700  13.609 4.02e-14 ***
## cyl         -1.58728    0.71184  -2.230   0.0337 *  
## disp        -0.02058    0.01026  -2.007   0.0542 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.055 on 29 degrees of freedom
## Multiple R-squared:  0.7596, Adjusted R-squared:  0.743 
## F-statistic: 45.81 on 2 and 29 DF,  p-value: 1.058e-09

On a déjà plus de détails.

On peut extraire divers éléments issus de lm() :

summary(ma_regression)$coef[,"Estimate"] # les beta chapeau

## (Intercept)         cyl        disp 
## 34.66099474 -1.58727681 -0.02058363

summary(ma_regression)$coef[,"Std. Error"] # les écarts-type

## (Intercept)         cyl        disp 
##  2.54700388  0.71184427  0.01025748

summary(ma_regression)$coef[,"Estimate"]["cyl"] # le beta chapeau de la variable "cyl"

##       cyl 
## -1.587277

summary(ma_regression)$coef[,"Estimate"][["cyl"]] # la même valeur, mais sans le nom associé, notez les doubles crochets

## [1] -1.587277

On va maintenant utiliser ces éléments pour construire une simulation de Monte-Carlo de la distribution de l’estimateur des MCO \(\hat\beta_x\) dans un cadre qui respecte les hypothèses de Gauss-Markov. On procède comme précédemment : une fonction qui simule un échantillon et retourne le \(\hat\beta_x\) estimé, puis replicate() pour répéter, et map() pour faire varier la taille d’échantillon.

library(tidyverse)

betareg<-function(n) {
  x <- rnorm(n,mean=0,sd=1)
  epsilon <- rnorm(n,mean=0,sd=1) # epsilon suit une loi normale non corrélée à x
  y <- 1+x+epsilon
  mareg <- lm(y~x)
  coeffs <- summary(mareg)$coef[,"Estimate"]
  betax <- coeffs[["x"]]
  betax # on retourne le beta chapeau de x
}

n_grid <- c(10, 100) # les tailles d'échantillon que l'on veut comparer

Comme il n’y a ici qu’un seul paramètre qui varie (n), on peut se contenter de map() (au lieu de pmap()). Pour chaque valeur de n, on lance replicate() afin d’obtenir 10000 estimations de \(\hat\beta_x\). On range le tout dans une data frame que l’on déplie avec unnest(), exactement comme dans la section précédente.

K <- 10000 # nombre de réplications par scénario

data_resultats <- tibble(n = n_grid) %>%
  mutate(betax = map(n, \(n) replicate(K, betareg(n)))) %>%
  unnest(betax)

head(data_resultats)

## # A tibble: 6 × 2
##       n betax
##   <dbl> <dbl>
## 1    10 0.395
## 2    10 1.06 
## 3    10 1.36 
## 4    10 0.813
## 5    10 0.355
## 6    10 1.52

data_resultats %>%
  group_by(n) %>%
  summarise(moy=mean(betax),ecty=sd(betax))

## # A tibble: 2 × 3
##       n   moy  ecty
##   <dbl> <dbl> <dbl>
## 1    10 1.000 0.386
## 2   100 1.00  0.100

ggplot(data=data_resultats,mapping=aes(x=betax,group=n,color=factor(n)))+geom_density()

Les résultats ci-dessus sont-ils ceux que l’on s’attendait à avoir ? Que nous dit la théorie ?

3.4 Pour résumer

La démarche générale d’une simulation de Monte-Carlo tient en quatre étapes :

Étape	Outil utilisé
Fonction de simulation	une fonction qui retourne directement la valeur
Grille des scénarios	`expand_grid(n = n_grid, ...)`
Lancer les réplications	`pmap(..., \(...) replicate(K, func(...)))`
Récupérer une `data frame` de résultats	`unnest()` sur la colonne-liste

L’avantage de cette approche est qu’elle ne repose que sur des fonctions de base de R et du tidyverse, que vous réutiliserez dans bien d’autres contextes.