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 :
library(boot)
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 :
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 :
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 :
names(repl.np)
dont les plus utiles sont :
repl.np$t
repl.np$t0
Pour utiliser l'approche bootstrap paramétrique, c'est légèrement plus compliqué :
Dans le cas de la moyenne empirique comme estimateur et des échantillons bootstrap tirés selon une loi normale, ça donne ça:
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
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$.
Pour le calcul de divers intervalles bootstrap, on utilise la fonction boot.ci :
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 :
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 :
names(icb.np)
Voici le niveau de confiance nominal $1-\alpha$ ainsi que les limites de l'intervalle $\mathcal I^*_{\text{norm}}$ :
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 :
icb.np$basic
Les mêmes informations pour l'intervalle $\mathcal I^*_{\text{perc}}$ :
icb.np$percent
et pour l'intervalle $\mathcal I^*_{\mathrm{BC_a}}$ :
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 :
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 :
boot.ci(repl.np.var)$stud
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.
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 ?