Le package boot

Le package boot contient des fonctions puissantes et optimisées pour le bootstrap. En effet, les fonctions que vous avez programmées au dernier notebook sont toutes disponibles dans le package boot. Cela ne veut pas dire que votre travail a été inutile. Au fait, la mise en oeuvre des différentes fonctions bootstrap vous a sûrement aidé a mieux comprendre les méthodes !

Maintenant chargons le package par l'instruction suivante :

In [ ]:
library(boot)

Répliques bootstrap

La fonction de base est la fonction boot. Elle génère des répliques bootstrap d'un estimateur $T$ - au choix par le bootstrap paramétrique ou non paramétrique. Elle a des nombreuses options :

boot(data, statistic, R, sim = "ordinary", stype = c("i", "f", "w"), strata = rep(1,n), L = NULL, m = 0, weights = NULL, ran.gen = function(d, p) d, mle = NULL, simple = FALSE, ..., parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), cl = NULL)

Voyons d'abord le bootstrap nonparamétrique. Les arguments de boot sont les suivants : le vecteur des données (data), une fonction qui définit l'estimateur $T$ (statistics) et le nombre $B$ de répliques souhaités. Par exemple :

In [ ]:
n <- 20
obs <- rnorm(n)
B <- 10
boot.mean <- function(obs,indice){
    moy.boot <- mean(obs[indice]) # permet à boot de choisir les données
    return(moy.boot)
}
repl.np <- boot(obs,boot.mean,B)
repl.np

On obtient toute une liste d'objets en sortie :

In [ ]:
names(repl.np)

dont les plus utiles sont :

  • repl\$t0 : la valeur $t=T(x_1,\dots,x_n)$ de l'estimateur sur les données observées
  • repl\$t : le vecteur des répliques bootstrap de l'estimateur
In [ ]:
repl.np$t
In [ ]:
repl.np$t0

Pour utiliser l'approche bootstrap paramétrique, c'est légèrement plus compliqué :

  • il faut préciser sim='parametric',
  • statistic est une fonction qui prend en argument un vecteur de données et renvoie la valeur de l'estimateur sur ces données,
  • ran.gen est la fonction qui génère un échantillon bootstrap paramétrique. Elle prend en argument les données et une liste de paramètres.
  • mle prend les valeurs observées des paramètres (que boot passera à la fonction ran.gen).

Dans le cas de la moyenne empirique comme estimateur et des échantillons bootstrap tirés selon une loi normale, ça donne ça:

In [ ]:
ran.gen.norm <- function(data,param){
  rnorm(length(data),param$mu,param$sig)
}
repl.param <- boot(obs,mean,B,sim='parametric',ran.gen=ran.gen.norm,mle=list(mu=mean(obs),sig=sd(obs)))
repl.param

Question 1

Considérons, une fois de plus, un échantillon $(X_1,\dots,X_n)$ de loi exponentielle $\mathcal E(\theta)$. La quantité à estimer est $1/\mathbb E[X_1]$, que l'on estime naturellement par $T=1/\bar X_n$.

  • Quelle est l'instruction pour créer 1000 répliques bootstrap non paramétrique de l'estimateur $T$ en utilisant la fonction boot ?
  • A l'image de ren.gen.norm, écrire une fonction ren.gen.exp pour la génération d'échantillons bootstrap de loi exponentielle.
  • Utiliser boot pour créer des répliques bootstrap paramétrique de $T$.

Intervalles de confiance

Pour le calcul de divers intervalles bootstrap, on utilise la fonction boot.ci :

boot.ci(boot.out, conf = 0.95, type = "all", index = 1:min(2,length(boot.out$t0)), var.t0 = NULL, var.t = NULL, t0 = NULL, t = NULL, L = NULL, h = function(t) t, hdot = function(t) rep(1,length(t)), hinv = function(t) t, ...)

Le premier argument est un objet obtenu par un appel à boot (essentiellement, ce sont des répliques bootstrap de l'estimateur).

Avec type on précise le(s) intervalles bootstrap souhaité(s) : au choix c("norm","basic", "stud", "perc", "bca") ou simplement type = "all" (par défaut) pour tous les intervalles.

Voici l'appel le plus simple de boot.ci (qui produit un warning car l'intervalle studentisé nécessite plus d'information) pour calculer des intervalles de confiance bootstrap non paramétrique de la moyenne de la loi normale basés sur la moyenne empirique :

In [ ]:
B <- 1000
repl.np <- boot(obs,boot.mean,B)
icb.np <- boot.ci(repl.np)
icb.np

Au fait, l'instruction boot.ci(repl.np) calcule les quatre intervalles bootstrap suivants : l'intervalle par approximation normale $\mathcal I^*_{\text{norm}}$, l'intervalle de base $\mathcal I^*_{\text{basic}}$, l'intervalle par la méthode de percentile de base $\mathcal I^*_{\text{perc}}$ et l'intervalle par la méthode des percentiles ajustée $\mathcal I^*_{\mathrm{BC_a}}$. Plus précisément, boot.ci renvoie une liste avec les éléments suivants :

In [ ]:
names(icb.np)

Voici le niveau de confiance nominal $1-\alpha$ ainsi que les limites de l'intervalle $\mathcal I^*_{\text{norm}}$ :

In [ ]:
icb.np$normal

Pour l'intervalle $\mathcal I^*_{\text{basic}}$ on a les informations suivantes : d'abord le niveau de confiance nominal, ensuite l'ordre utilisé des statistiques d'ordres dans les limites de l'intervalle, et enfin les limites de l'intervalle bootstrap :

In [ ]:
icb.np$basic

Les mêmes informations pour l'intervalle $\mathcal I^*_{\text{perc}}$ :

In [ ]:
icb.np$percent

et pour l'intervalle $\mathcal I^*_{\mathrm{BC_a}}$ :

In [ ]:
icb.np$bca

Pour le calcul de l'intervalle bootstrap studentisé $\mathcal I^*_{\text{stud}}$ il faut fournir des estimations de la variance de $T$ pour les différents échantillons bootstrap. Une façon de faire consiste à modifier la fonction boot.mean pour qu'elle renvoie également l'estimateur de la variance :

In [ ]:
obs <- rnorm(20,10,2)       
boot.mean.var <- function(obs,indice){
    obs.boot <- obs[indice] 
    moy.boot <- mean(obs.boot) 
    var.boot <- var(replicate(50,mean(sample(obs.boot,replace=T))))
    return(c(moy.boot,var.boot))
}
repl.np.var <- boot(obs,boot.mean.var,R=1000)
repl.np.var

Maintenat on obtient l'intervalle bootstrap studentisé par l'instruction suivante :

In [ ]:
boot.ci(repl.np.var)$stud

Question 2

  • Pour un jeu de données simulé de loi exponentielle $\mathcal E(\theta)$ calculer les intervalles de confiance bootstrap $\mathcal I^*_{\text{norm}}$, $\mathcal I^*_{\text{basic}}$, $\mathcal I^*_{\text{perc}}$ et $\mathcal I^*_{\mathrm{BC_a}}$ pour $1/\mathbb E[X_1]$ en utilisant des répliques de $T=1/\bar X_n$ par le bootstrap non paramétrique.
  • Mettre en oeuvre le calcul de l'intervalle studentisés par le bootstrap non paramétrique.
  • Quant au bootstrap paramétrique calculer les intervalles de confiance bootstrap $\mathcal I^*_{\text{norm}}$, $\mathcal I^*_{\text{basic}}$ et $\mathcal I^*_{\text{perc}}$.
  • Observer que boot.ci ne fournit pas l'intervalle $\mathcal I^*_{\mathrm{BC_a}}$. Au fait, boot.ci ne sait le calculer que dans le cas non paramétrique. Dans le cas paramétrique, la définition de l'intervalle est différente, mais elle n'est pas implementée dans boot.ci. Nous allons nous contenter des autres intervalles.
  • Mettre en oeuvre le calcul de l'intervalle studentisé par le bootstrap paramétrique.

Question 3

Ecrire un programme pour simuler un grand nombre de jeu de données de loi exponentielle et calculer tous les intervalles de confiance $\mathcal I^*_{\text{norm}}$, $\mathcal I^*_{\text{basic}}$, $\mathcal I^*_{\text{stud}}$ et $\mathcal I^*_{\text{perc}}$ pour $1/\mathbb E[X]$ par les deux méthodes : bootstrap paramétrique et non paramétrique. Enfin estimer le niveau de confiance effectif de chaque intervalle. Commenter les résultats.

Question 4

On veut étudier ce qui se passe si les données ne suivent pas exactement une loi exponentielle. Maintenant les données suivent la loi Gamma $\Gamma(\alpha,\theta)$. La quantité inconnue est toujours $1/\mathbb E[X]=\theta/\alpha$, estimée par $T=1/\bar X_n$.

Modifier votre programme pour générer des données de loi Gamma. Ensuite, appliquer exactement les mêmes calculs comme si les données avaient une loi exponentielle (notamment l'échantillons bootstrap paramétrique est toujours tiré selon une loi exponentielle).

Rappelons que $\mathcal E(\theta)=\Gamma(1,\theta)$. Comment évolue les différent niveaux de confiance effectif lorsque les données suivent une loi $\Gamma(\alpha,\theta)$ avec des valeurs de $\alpha$ inférieures à 1 ?