POL6022 - Méthodes quantitatives avancées

Bootstrapping

En statistique, on veut souvent rééchantillonner pour tester la précision de nos estimations d’échantillon. On appelle cela le bootstrapping — une méthode fondée sur de multiples tirages aléatoires avec remise. Dans ce labo, nous allons : (1) générer des échantillons aléatoires à partir de fonctions de distribution (avec remise), (2) introduire des boucles et (3) calculer l’écart-type de la distribution de nos échantillons bootstrap.

Dans ce matériel, je m’appuie sur le code mis à disposition par Carsey et Harden (2014 : 217).

Fonctions pertinentes : rnorm(), rbernoulli(), rbinom(), rpois(), for()

1. Génération de données

1.1 Distribution normale

La fonction rnorm() permet de générer une suite aléatoire de nombres à partir d’une distribution normale hypothétique centrée sur une moyenne de 20 et un écart-type de 4,5. Ici, nous choisissons de générer 2000 observations.

set.seed(300) # Définition du "seed" pour la reproductibilité

myData <- rnorm(n=2000, mean=20, sd=4.5) # Création d'une distribution normale aléatoire

Pour confirmer que le vecteur myData a été créé correctement, on vérifie que (1) il y a 2000 observations (length), (2) la moyenne est approximativement 20 (mean) et (3) l’écart-type est approximativement 4,5 (sd). Note : on ne s’attend pas à obtenir exactement les valeurs entrées dans rnorm, puisque la fonction introduit de l’aléatoire.

length(myData) # Combien d’observations ?

## [1] 2000

mean(myData) # Quelle est la moyenne ?

## [1] 20.25773

sd(myData) # Quel est l’écart-type ?

## [1] 4.590852

À titre visuel, voici un graphique du vecteur myData. La ligne rouge pointillée représente la moyenne de cette distribution.

1.2 Autres distributions

Notez que vous pouvez créer une variété de distributions dans R. Voici comment créer celles dont nous avons discuté en classe.

# Définition du "seed" pour la reproductibilité
set.seed(300)

# Chargement du package purrr
library(purrr)

# Création d'une variable de Bernoulli aléatoire avec n=1 (1 tirage)
# et p=0.5 (probabilité de succès de 50 %)
rbernoulli(n=1, p=0.5)

## Warning: `rbernoulli()` was deprecated in purrr 1.0.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## [1] TRUE

# N'oubliez pas de définir votre "seed" à 123!
set.seed(123)

# Vous devez entrer les bons arguments dans cette fonction :
# rbinom(n=?, size=?, prob=?)

# Réponse

set.seed(123)

rbinom(n=50, size=10, prob=0.5)

##  [1] 4 6 5 7 7 2 5 7 5 5 8 5 6 5 3 7 4 2 4 8 7 6 6 9 6 6 5 5 4 3 8 7 6 6 2 5 6 4
## [39] 4 4 3 5 5 4 3 3 4 5 4 7

set.seed(300) # Définition du "seed" pour la reproductibilité

rpois(n=14, lambda=30) # Création d'une variable de Poisson aléatoire (n=14, lambda=30)

##  [1] 37 34 32 33 29 38 34 32 36 31 42 29 22 36

2. Boucles

Introduisons les boucles avec un exemple très simple. Ici, on veut afficher 5 observations tirées de myData à l’aide de sample() et print(). Autrement dit, on tire et on affiche une observation, cinq fois de suite.

set.seed(300) # Définition du "seed" pour la reproductibilité

for (i in 1:5) # Spécification du nombre d’itérations
{
  obs <- sample(myData, size=1) # Tirage d'une observation de myData
  print(obs) # Affichage de cette observation
  cat("J’ai terminé", i, "itérations \n") # Message après chaque itération
}

## [1] 13.85058
## J’ai terminé 1 itérations 
## [1] 18.4574
## J’ai terminé 2 itérations 
## [1] 18.57355
## J’ai terminé 3 itérations 
## [1] 22.5185
## J’ai terminé 4 itérations 
## [1] 22.79222
## J’ai terminé 5 itérations

On peut utiliser des boucles pour automatiser un grand nombre de choses. Dans les prochaines sous-sections, nous les utiliserons pour illustrer deux principes importants : la loi des grands nombres et le théorème central limite.

2.1 Loi des grands nombres

Créons une fonction die_mean() pour calculer la moyenne de n lancers de dé.

die <- 1:6 # Création d'un vecteur "die" avec tous les nombres de 1 à 6 

die_mean <- function(n) {
  # Échantillonnage avec remise "n" fois dans le vecteur die
  mean(sample(die, size=n, replace=TRUE))
}

Calcul de la moyenne de n lancers de dé avec die_mean().

set.seed(500) # Définition du "seed" pour la reproductibilité

for (n in c(10,100,1000,10000,100000))
{
  result <- die_mean(n)
  cat("La moyenne de", n, "lancers de dé est", result, "\n")
}

## La moyenne de 10 lancers de dé est 3.2 
## La moyenne de 100 lancers de dé est 3.44 
## La moyenne de 1000 lancers de dé est 3.492 
## La moyenne de 10000 lancers de dé est 3.4994 
## La moyenne de 1e+05 lancers de dé est 3.50887

La moyenne d’échantillon se rapproche de la valeur attendue (paramètre de population) à mesure que la taille de l’échantillon augmente. C’est la loi des grands nombres.

2.2 Théorème central limite

Tirer i observations à partir d’un vecteur de moyennes de 15 lancers de dé avec replicate() et tracer les résultats avec hist().

set.seed(100) # Définition du "seed" pour la reproductibilité

for (i in c(100,1000,100000))
{
  x <- replicate(i, {
    mm <- die_mean(15) # Moyenne de 15 lancers
    mean(mm)
  })
  
  a <- round(sd(x), digits=3)
  
  hist(x,
       xlab="Moyenne de 15 lancers de dé",
       xlim=c(1,6),
       col="#4286f4",
       main=paste("Histogramme pour ", i, " tirages aléatoires (sd = ", a, ")", sep=""))
}

Vous pouvez cliquer sur ces graphiques pour les examiner de plus près !

La distribution de la moyenne d’échantillon tend vers une distribution normale lorsque la taille de l’échantillon augmente. C’est le théorème central limite.

3. Bootstrapping

Le bootstrapping est une façon de tenir compte de l’incertitude lorsque nous mesurons des estimations d’échantillon. En utilisant des tirages aléatoires avec remise, on calcule nos estimations sur de nombreux échantillons aléatoires tirés de notre distribution d’origine.

Ici, nous calculons la moyenne de 1000 échantillons tirés d’une distribution normale aléatoire. La moyenne de chacun de ces échantillons est stockée dans un vecteur (bootstrap.results). Ce qui nous intéresse, c’est l’incertitude de l’estimation (ici : la moyenne), mesurée par l’écart-type de la distribution des moyennes.

set.seed(200) # Définition du "seed" pour la reproductibilité

n.samples <- 1000 # Nombre d’échantillons bootstrap

bootstrap.results <- c() # Création d'un vecteur vide pour stocker les résultats

for (i in 1:n.samples)
{
  bootstrap.results[i] <- mean(rnorm(2000, 20, 4.5)) # Moyenne de l’échantillon bootstrap
}

length(bootstrap.results) # Vérification : 1000 moyennes

## [1] 1000

summary(bootstrap.results) # Vérification

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.64   19.93   20.00   20.00   20.07   20.32

sd(bootstrap.results) # Écart-type de la distribution des moyennes

## [1] 0.1041927

par(mfrow=c(2,1), pin=c(5.8,0.98)) # Combinaison de graphiques

hist(bootstrap.results,
     col="#d83737",
     xlab="Moyenne",
     main=paste("Moyennes de 1000 échantillons bootstrap (DGP)"))

hist(myData,
     col="#37aad8",
     xlab="Valeur",
     main=paste("Distribution de myData"))

# N'oubliez pas de définir votre "seed" à 120!
set.seed(120)

# Écrivez votre code ici

# Réponse

set.seed(120) # Définition du "seed" pour la reproductibilité

n.samples <- 50 # Nombre d’échantillons bootstrap

bootstrap.results <- c() # Création d'un vecteur vide pour stocker les résultats

for (i in 1:n.samples)
{
  bootstrap.results[i] <- mean(rpois(n=200, lambda=7)) # Moyenne de l’échantillon bootstrap
}

summary(bootstrap.results) # Résultat

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   6.610   6.860   6.987   6.999   7.120   7.455