Dans ce TP on considère des données concernant un système de climatisation. Les observations sont les durées (en heure) entre deux pannes de la climatisation (source: Bootstrap methods and their application, Davison & Hinkley, 1997). On s'intèresse à l'estimation du taux de défaillance. Celui est estimé par $T=1/\bar X_n$ où $X_i$ désignent les observations. Voici les données :
aircond <- c(487,18,85,91,5,130,230,43,98,7,100,3)
Tracer l'histogramme des données aircond. Supposons que les données suivent une loi exponentielle $\mathcal E(\theta)$. Estimer le paramètre $\theta$ à partir des données et superposer la densité estimée à l'histogramme. Le choix de la loi exponentielle vous semble adéquat?
Notons $x_1,\dots,x_n$ des réalisations i.i.d. d'une loi $F$.
Rappelons que l'on obtient un échantillon bootstrap non paramétrique en tirant avec remise des valeurs dans l'échantillon $(x_1,\dots,x_n)$.
Si on a un modèle paramétrique des données en sorte que $X_i\sim F_\theta$ avec un paramètre $\theta$ inconnu et si on dispose d'un estimateur $T$ de $\theta$, on créé un échantillon bootstrap paramétrique en simulant de la loi estimée $F_t$ où $t=T(x_1,\dots,x_n)$ est la valeur de l'estimateur sur les données observées.
Selon le type d'échantillon bootstrap utilisé (paramétrique ou non), on obtient des répliques bootstrap paramétrique ou non paramétrique $T^*$ de l'estimateur $T$. Autrement dit, $T^*$ est la valeur de l'estimateur $T$ sur un échantillon bootstrap.
boot.np <- function(obs,R){
return(replic.boot)
}
boot.exp <- function(obs,R){
return(replic.boot)
}
estim.exp <- function(theta,n,R){
return(estim)
}
Soit $X_1,\dots,X_n$ un échantillon de variables i.i.d. de loi uniforme $U(\theta, \theta+ 1)$ pour $\theta\in\mathbb R$.
On considère l'estimateur $T_n = X_{(1)} = \min\{X_i,\dots, X_n\}$ de $\theta$. On s'intéresse à la distribution de la statistique $$L_n = n(T_n − \theta).$$
On peut prouver rigoureusement le résultat asymptotique suivant : $$L_n\stackrel{\mathcal L}{\longrightarrow}\mathcal E(1),\qquad\text{lorsque }n\to\infty.$$
Vérifier à l'aide des simulations à partir de quelle taille d'échantillon $n$ l'approximation limite donnée ci-dessus est bonne. Pour cela :
simul_min <- function(theta,n,R){
}
Maintenant nous voulons faire du bootstrap non paramétrique afin d'approcher la loi de la statistique $L_n$.
Ecrire une fonction boot_min similaire à la fontion simul_min qui prend en argument un échantillon et le nombre $R$ d'échantillons bootstrap à générer. Cette fonction trace l'histogramme et le QQ-plot des répliques bootstrap $L_n^*$. Qu'observez-vous en appelant cette fonction ? Obtient-on les mêmes résultats que par des simuations de la question précédente ? Interprétez.
boot_min <- function(obs,R){
}
Pour mieux comprendre ce qui se passe en effectuant du bootstrap non paramétrique, vérifier le nombre d’éléments différents dans le vecteur de répliques bootstrap $L^*_{n,r}, r=1,\dots,R$? Regardez pour des grandes valeurs de $n$ et de $R$. Expliquez.
Est-ce que le bootstrap paramétrique fait mieux ?
Mettre en oeuvre le bootstrap paramétrique où l'échantillon bootstrap est tiré selon la loi uniforme $U[t,t+1]$ où $t=T(x_1,\dots,x_n)$ est la valeur observée de l'esimtateur $T$ sur les données. Tracer l'histogramme et le QQ-plot correspondant. Interpréter.
Dans cette partie nous considérons des réalisations $x_1,\dots,x_n$ i.i.d. de loi exponentielle $\mathcal E(\theta)$ avec $\theta>0$, et l'estimateur $T=1/\bar X_n$. Notons $t=T(x_1,\dots,x_n)$ la valeur observée de l'estimateur sur les données.
Nous allons mettre en oeuvre les différents intervalles de confiance par le bootstrap (ICB) vus en cours.
Rappelons que le niveau de confiance nominatif d'un ICB est de $1-\alpha$ pour $\alpha$ fixé. Nous allons effectuer des simulations de Monte Carlo afin d'estimer le niveau de confiance effectif d'un ICB $\mathcal I^*$, c'est-à-dire déterminer la quantité $$\mathbb P_\theta(\theta\in\mathcal I^*)$$ et comparer cette valeur à la valeur nominative de $1-\alpha$ afin de juger la qualité de l'intervalle.
L'ICB de niveau $1-\alpha$ qui repose sur une approximation par la loi normale est donné par $$\mathcal I^*_{\text{norm}} = \left[t-\beta_B^*-\sqrt{v^*_B}z_{1-\alpha/2},t-\beta_B^*-\sqrt{v^*_B}z_{\alpha/2} \right],$$ où $z_{\alpha/2}$ et $z_{1-\alpha/2}$ désignent les quantiles de la loi normale standard $\mathcal N(0,1)$ et $\beta_B^*$ et $v^*_B$ sont les estimateurs bootstrap du biais et de la variance de l'estimateur $T$ définis par $$\beta_B^*=\bar T^*-t,\qquad v^*_B = \frac1B\sum_{b=1}^B(T^{*(b)}-\bar T^*)^2,\qquad \text{ où } \quad\bar T^*=\frac1B\sum_{b=1}^BT^{*(b)},$$ et $T^{*(b)},b=1,\dots,B$ sont $B$ répliques bootstrap (paramétrique ou non paramétrique) de $T$.
biais.var.boot <- function(obs,B,param=FALSE){
return(list(biais = b, var = v))
}
icb.norm <- function(obs,alpha=.05,B=300,param=FALSE){
return(icb)
}
niveau.icbnorm <- function(theta,n,M,alpha=0.05,param=FALSE){
return(prop)
}
L'intervalle de confiance par bootstrap de base est définit par $$\mathcal I^*_{\text{basic}} = \left[2t-T^*_{(\gamma_1)}, 2t-T^*_{(\gamma_2)} \right],$$ où $T^*_{(s)}$ désigne la $s$-ième statistique d'ordre associée aux répliques bootstrap $T^*_1,\dots,T^*_R$ et $$\gamma_1=\left\lceil\left(1-\frac\alpha 2\right)R\right\rceil\qquad\text{ et }\qquad\gamma_2=\left\lceil \frac\alpha2R\right\rceil.\qquad\qquad\qquad(*)$$
L'ICB studentisé est définit par $$\mathcal I^*_{\text{stud}} = \left[t-\sqrt{v_R^*}Z^*_{(\gamma_1)}, t-\sqrt{v_R^*}Z^*_{(\gamma_2)} \right],$$ avec $\gamma_1$ et $\gamma_2$ donnés en $(\ast)$ et $$Z^*_r = \frac{T^*_r-t}{\sqrt{v^*_{r,B}}},$$ où $v^*_{r,B}$ désigne l'estimateur bootstrap de la variance de $T$ associée au $r$-ième échantillon bootstrap $x^*_{r,1},\dots,x^*_{r,n}$ (double bootstrap).
L'ICB par la méthode des percentiles de base est définit par $$\mathcal I^*_{\text{perc}} = \left[T^*_{(\gamma_2)}, T^*_{(\gamma_1)} \right].$$
Au vu de vos résultats de simulation comparer la performance des différents types d'intervalles de confiance et expliquer leur différences. Quelle approche, paramétrique ou non paramétrique, donne des meilleurs performances (dans le cadre étudié ici) ?
Calculer les différents intervalles de confiance pour les données aircond et comparer-les.