TP1 Méthode de Bootstrap¶

1) Motivation¶

Considérons qu'on dispose d'un échantillon $(X_1,\dots,X_n)$ d'une certaine loi $\mathbb{P}_{\theta^\ast}$ qui dépend d'un paramètre $\theta^\ast$ inconnu. Au travers d'exemples, nous allons voir comment estimer $\theta^\ast$ à partir des observations et comment évaluer la qualité et la précision d'un estimateur.

On prendra comme exemple la loi normale et la loi Gamma. Ces deux lois dépendent de deux paramètres:

  • pour la loi normale $\theta = (\mu,\sigma)$ dont la densité est $$ f(x;\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-(x-\mu)^2/(2\sigma^2)}$$
  • pour la loi Gamma $\theta = (\alpha,\beta)$ dont la densité est $$ f(x;k,\beta) = \frac{1}{\Gamma(\alpha)\beta^k}x^{k-1}e^{-x/\beta }\mathbf{1}_{]0,+\infty[}(x).$$

2) Le cas d'un échantillon Gaussien¶

  • Si $X\sim \mathcal{N}(\mu,\sigma^2)$, on a $\mathbb{E}(X)=\mu$ et $var(X) = \sigma^2$, ce qui nous ammène aux estimateurs

$$ \bar X_n = \frac{1}{n}\sum_{j=1}^n X_j \qquad \text{et}\qquad \bar \sigma^2_n = \frac{1}{n-1}\sum_{j=1}^n (X_j - \bar X_n)^2,$$

de respectivement $\mu$ et $\sigma^2$. La consistance de ces estimateurs est garantie par la loi des grands nombres: presque sûrement

$$\lim_{n\to +\infty} \bar X_n = \mu \qquad\text{et}\qquad \lim_{n\to +\infty} \sigma^2_n = \sigma^2.$$

Concernant la qualité de l'estimation, on remarque tout d'abord que $\bar X_n \sim \mathcal{N}(\mu,\sigma^2/n)$, ou autrement dit

$$ N = \frac{\sqrt{n}}{\sigma}(\bar X_n - \mu) \sim \mathcal{N}(0,1).$$

De plus,

$$ C = \frac{n-1}{\sigma^2} \bar\sigma_n^2 \sim \chi^2_{n-1}$$

où $\chi^2_{n-1}$ représente la loi du $\chi^2$ à $n-1$ degré de liberté. Ainsi,

$$ \frac{\sqrt{n}}{\bar \sigma_n}(\bar X_n - \mu) = \frac{N}{\sqrt{C/(n-1)}}\sim t(n-1) $$

où $t(n-1)$ représente la loi de Students à $n-1$ degrés de liberté.

a) Pour $\mu$.¶

Pour ce paramètre, en notant $q_{n-1,l}$ le quantile d'ordre $l$ de la loi $t(n-1)$ (autrement dit $F_{t(n-1)}(q_{n-1,l}) = l$), on remarque que

$$ \mathbb{P}\Big(\mu \in \Big[\bar X_n - q_{n-1,1-\alpha/2}\frac{\bar \sigma_n}{\sqrt{n}},\bar X_n + q_{n-1,1-\alpha/2}\frac{\bar \sigma_n}{\sqrt{n}}\Big]\Big) = \mathbb{P}\Big( \Big|\frac{\sqrt{n}}{\bar \sigma_n}(\bar X_n - \mu)\Big|\leq q_{n-1,1-\alpha/2}\Big) = \mathbb{P}(|T|\leq q_{n-1,1-\alpha/2}) $$

où $T\sim t(n-1)$ et $\alpha\in]0,1[$. Or, on a

$$\mathbb{P}(|T|\leq q_{n-1,1-\alpha/2}) = 2 F_{t(n-1)}(q_{n-1,1-\alpha/2}) - 1 =1-\alpha.$$

En conclusion,

$$ \mathbb{P}\Big(\mu \in \Big[\bar X_n - q_{n-1,1-\alpha/2}\frac{\bar \sigma_n}{\sqrt{n}},\bar X_n + q_{n-1,1-\alpha/2}\frac{\bar \sigma_n}{\sqrt{n}}\Big]\Big) = 1-\alpha $$

On parle d'intervalle de confiance de niveau $1-\alpha$, et $\alpha$ représente la probabilité que $\mu$ ne soit pas dans l'intervalle. Les valeurs typiques de $\alpha$ sont $0.05$ ou $0.01$.

Regardons un exemple numérique, avec $\alpha=0.05$, $\mu=0$, $\sigma^2=2$ et un échantillon de taille $n=50$, et regardons le comportement des estimateurs et de l'intervalle de confiance.

In [1]:
import numpy as np
import scipy.stats as sps 
import matplotlib.pyplot as plt
from IPython.core.pylabtools import figsize
In [2]:
n = 50
mu = 0.
sig = np.sqrt(2)
alpha = 0.05

X = sps.norm(mu, sig).rvs(n)

On peut calculer l'estimation de $\mu$

In [3]:
mu_hat = X.mean()
mu_hat
Out[3]:
0.2701274214575453

et de l'écart type

In [4]:
sig_hat = X.std()
sig_hat
Out[4]:
1.5049667537317222

et calculer un intervalle de confience

In [5]:
q = sps.t(n-1).ppf(1-alpha/2) #la méthode ppf correspond au pourcentile/quantile

IC1_mu = mu_hat - q * sig_hat / np.sqrt(n)
IC2_mu = mu_hat + q * sig_hat / np.sqrt(n)

print("les bornes de l'IC sont %s et %s" % (IC1_mu, IC2_mu))

les bornes de l'IC sont -0.15757939757862666 et 0.6978342404937172

On peut maintenant répéter ce calcul pour voir jouer le rôle de $\alpha$. Ici $\alpha=0.05$. On peut donc s'attendre à ce qu'en moyenne la vraie valeur de $\mu$ ne soit pas dans l'intervalle de confiance 1 fois sur 20. On va faire 40 répétitions et tracer les intervalles de confiance.

In [6]:
def plot_IC(mu, sig, n, alpha, M):
    
    """représente graphiquement M IC de niveaux 1- alpha
    pour mu à partir d'un échantillon de taille n"""
    
    q = sps.t(n-1).ppf(1-alpha/2)
    
    for j in range(M):
        
        X = sps.norm(mu, sig).rvs(n)
        mu_hat = X.mean()
        sig_hat = X.std()
        IC1 = mu_hat - q * sig_hat / np.sqrt(n)
        IC2 = mu_hat + q * sig_hat / np.sqrt(n)
        
        plt.vlines(j, IC1, IC2)
        plt.hlines(mu, -1, M)
In [7]:
M = 40
plot_IC(mu, sig, n, alpha, M)
figsize(2,5)
No description has been provided for this image

b) Pour $\sigma^2$¶

Pour ce paramètre, en notant encore $q_{n-1,l}$ le quantile d'ordre $l$ de cette fois de la loi $\chi^2(n-1)$ , on remarque que

$$ \mathbb{P}\Big(\sigma^2 \in \Big[\frac{n-1}{q_{n-1,1-\alpha/2}}\bar \sigma_n^2, \frac{n-1}{q_{n-1,\alpha/2}}\bar \sigma_n^2\Big]\Big) = \mathbb{P}\Big(\frac{n-1}{\sigma^2} \bar\sigma^2_n \in [q_{n-1,\alpha/2}, q_{n-1,1-\alpha/2}]\Big) = \mathbb{P}\Big(C \in [q_{n-1,\alpha/2}, q_{n-1,1-\alpha/2}]\Big) $$

où $C\sim \chi^2(n-1)$ et $\alpha\in]0,1[$. Or, on a

$$\mathbb{P}\Big(C \in [q_{n-1,\alpha/2}, q_{n-1,1-\alpha/2}]\Big)= F_{\chi^2(n-1)}(q_{n-1,1-\alpha/2}) - F_{\chi^2(n-1)}(q_{n-1,\alpha/2}) =1-\alpha/2 - \alpha/2 = 1-\alpha.$$

On dit alors que l'intervalle

$$ \Big[\frac{n-1}{q_{n-1,1-\alpha/2}}\bar \sigma_n^2, \frac{n-1}{q_{n-1,\alpha/2}}\bar \sigma_n^2\Big] $$

de confiance de niveau $1-\alpha$ pour $\sigma^2$. Déterminons le pour notre échantillon.

In [8]:
q1 = sps.chi2(n-1).ppf(alpha/2)
q2 = sps.chi2(n-1).ppf(1-alpha/2)

IC1_sig2 = (n - 1) * sig_hat**2 / q2 
IC2_sig2 = (n - 1) * sig_hat**2 / q1

print("les bornes de l'IC sont %s et %s" % (IC1_sig2, IC2_sig2))

les bornes de l'IC sont 1.5804259057125298 et 3.5170849427966315

Dans le cas d'un échantillon Gaussien, pour évaluer la précision de l'estimation via les intervalles de confiance, on peut faire appel directement à des lois explicites. Ici pour $\mu$ et $\sigma^2$, on arrive à obtenir des quantités indépendantes de l'échantillon et des paramètres qui nous permettent d'utiliser les quantiles de la loi de Students ou du $\chi^2$.

Dans le cas général, il n'est pas forcement possible de faire de telles manipulations. Comment faire alors? Une porte de sortie est l'utilisation du théorème central limite comme on le verra plus bas avec le cas d'un échantillon de loi Gamma. Une autre possibilité est d'utiliser la méthode du Bootstrap. Nous allons voir deux approche de cette méthode: Le bootstrap paramétrique et non-paramétrique.

La technique du bootstrap date des années 70 et a été inventée par Bradley Efron. Cette technique permet d'analyser la sensibilité des estimations en fonction du jeu de données, elles est uniquement basée sur la réplication du jeu de données original (rééchantillonnage).

3) Le Bootstrap (paramétrique)¶

Supposons qu'on dispose d'un échantillon $\mathbf{X} = (X_1,\dots,X_n)$ d'une certaine loi $\mathbb{P}_{\theta^\ast}$ qui dépend d'un paramètre $\theta^\ast$ inconnu, et supposons qu'on ait une estmation $\hat \theta$ de ce paramètre.

L'idée de la méthode consiste à

  • créer de nouveaux échantillons $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(B)}$ (bootstrap sample) à partir de la loi $\mathbb{P}_{\hat \theta}$
  • puis à partir de ces échantillons procéder à de nouvelles estimations $\hat \theta^{(1)},\dots,\hat \theta^{(B)}$ (bootstrap replicates)

a) Distribution Bootstrap¶

Avec l'échantillons $\hat \theta^{(1)},\dots,\hat \theta^{(B)}$ on peut construire la loi bootstrap de l'estimateur $\hat \theta$ à l'aide de la fonction de répartition empirique

$$ F_B(x) = \frac{1}{B}\sum_{j=1}^B \mathbf{1}_{\big(\hat \theta^{(j)}\leq x\big)}. $$

b) Intervalle de confiance Bootstrap¶

On peut aussi construire des intervalles de confiance de niveau $1-\alpha$. Pour cela, on ordonne de manière croissante l'échantillon $\hat \theta^{(1)},\dots,\hat \theta^{(B)}$ pour obtenir une famille vérifiant $\tilde \theta^{(1)}<\dots<\tilde \theta^{(B)}$. On retire la proportion $\alpha$ des plus petits échantillons et la proportion $1-\alpha/2$ des plus grands. Autrement dit les éléments $$ \tilde \theta^{([\alpha B/2])} \quad\text{et}\quad \tilde \theta^{([(1-\alpha/2) B])} $$ représentent les bornes de l'intervalle. Ici $[\cdot]$ correspond à la partie entière.

c) Estimation Bootstrap du biais et de la variance¶

Ces estimations sont respectivement données par

$$ \widehat{Biais} =\bar{\hat \theta} - \hat \theta\qquad\text{et}\qquad \widehat{Var} = \frac{1}{B-1}\sum_{j=1}^B \big(\hat\theta^{(j)} - \bar{\hat \theta}\big)^2 $$

avec

$$ \bar{\hat \theta} = \frac{1}{B}\sum_{j=1}^B \hat \theta^{(j)}.$$

A noter que dans cette méthode il y a deux sources d'erreur. Une due au fait qu'on utilise $\hat \theta$ et non $\theta^\ast$ (qu'on ne connait pas) pour rééchantillonner. La deuxième vient du fait qu'on utilise un nombre fini $B$ d'échantillon. On a donc là une erreur de type Monte-Carlo.

d) Revenons à notre exemple avec la loi normale.¶

Commençons par génrer la distribution Bootstrap de $\hat \theta$ de taille $B=20000$.

In [9]:
def dist_boot(data, B):
    
    """génèe un échantillon bootstrap pour mu et sigma de taille B"""
    
    mu_hat = data.mean()
    sig_hat = data.std()
    
    n = len(data)
    
    X = sps.norm(mu_hat, sig_hat).rvs((B, n))
    
    mu_boot = X.mean(axis=1)
    sig_boot = X.std(axis=1)
    
    return mu_boot, sig_boot
In [10]:
B = 20000

mu_boot, sig_boot = dist_boot(X, B)

mu_tilde = np.sort(mu_boot)
sig2_tilde = np.sort(sig_boot**2)

Représentons graphiquement l'histogramme normalisée et la fonction de répartition.

In [11]:
figsize(9, 9)

plt.subplot(221)

plt.hist(mu_boot, bins=80, density=True, alpha=0.2, label="échantillon");
plt.vlines(mu, 0, 2.5, label="vraie valeur")
plt.title("histogramme échantillon bootstrap pour $\mu$")
plt.legend()

plt.subplot(222)

y = np.linspace(1/B, 1, B)
plt.plot(mu_tilde, y)
plt.title("fonction de répartition bootstrap pour $\mu$")

plt.subplot(223)

plt.hist(sig_boot**2, bins=80, density=True, alpha=0.2, label="échantillon");
plt.vlines(sig**2, 0, 1, label="vraie valeur")
plt.title("histogramme échantillon bootstrap pour $\sigma^2$")
plt.legend()

plt.subplot(224)

y = np.linspace(1/B, 1, B)
plt.plot(sig2_tilde, y)
plt.title("fonction de répartition bootstrap pour $\sigma^2$");
No description has been provided for this image

On peut calculer les bornes pour les IC Bootstrap. On va considérer encore $\alpha = 0.05$.

In [12]:
IC1_mu_boot = mu_tilde[int(np.floor(alpha*B/2))]
IC2_mu_boot = mu_tilde[int(np.floor((1-alpha/2)*B))]

print("les bornes de l'IC pour mu sont %s et %s" % (IC1_mu_boot, IC2_mu_boot))
print("alors que celles calculées explicitement sont %s et %s" % (IC1_mu, IC2_mu))

les bornes de l'IC pour mu sont -0.15195829911475456 et 0.6909207740983553
alors que celles calculées explicitement sont -0.15757939757862666 et 0.6978342404937172
In [13]:
IC1_sig2_boot = sig2_tilde[int(np.floor(alpha*B/2))]
IC2_sig2_boot = sig2_tilde[int(np.floor((1-alpha/2)*B))]

print("les bornes de l'IC pour sigma^2 sont %s et %s" % (IC1_sig2_boot, IC2_sig2_boot))
print("alors que celles calculées explicitement sont %s et %s" % (IC1_sig2, IC2_sig2))

les bornes de l'IC pour sigma^2 sont 1.4275222325481107 et 3.189112308195373
alors que celles calculées explicitement sont 1.5804259057125298 et 3.5170849427966315

On voit que les résultats sont assez similaires sur cette exemple. On peut aussi estimer le biais et la variance Bootstrap

In [14]:
mu_boot_hat = mu_boot.mean()

biais_mu_boot = mu_boot_hat - mu_hat
var_mu_boot = np.mean( (mu_boot - mu_boot_hat)**2 )

print("le biais Bootstrap pour mu est %s" % biais_mu_boot)
print("la variance Bootstrap pour mu est %s" % var_mu_boot)
print("l'écart type pour mu est %s" % np.sqrt(var_mu_boot))

le biais Bootstrap pour mu est 0.00031702053617094217
la variance Bootstrap pour mu est 0.04543845472729605
l'écart type pour mu est 0.21316297691507324
In [15]:
sig2_boot_hat = np.mean(sig_boot**2)

biais_sig2_boot = sig2_boot_hat - X.var()
var_var_boot = np.mean( (sig_boot**2 - sig2_boot_hat)**2 )

print("le biais Bootstrap pour sigma^2 est %s" % biais_sig2_boot)
print("la variance Bootstrap pour sigma^2 est %s" % var_var_boot)
print("l'écart type pour sigma^2 est %s" % np.sqrt(var_var_boot))

le biais Bootstrap pour sigma^2 est -0.04075432194266604
la variance Bootstrap pour sigma^2 est 0.20346919464805688
l'écart type pour sigma^2 est 0.45107559748678144

4) Le cas d'un échantillon de loi Gamma¶

Prennons maintenant un exemple pour lequelle des calculs explicites ne sont pas directement accessibles pour évaluer la précision des estimateurs. On va considérer un échantillon de loi Gamma $\theta^\ast = (\alpha^\ast,\beta^\ast)$ dont la densité est $$ f(x;k,\beta) = \frac{1}{\Gamma(k)\beta^k}x^{k-1}e^{- x/\beta}\mathbf{1}_{]0,+\infty[}(x)$$

In [16]:
n = 2000
k = 1.5
beta = 2

X = sps.gamma(k, scale=beta).rvs(n)

Pour estimer $k$ et $\beta$, on va utiliser les relations

$$ \mathbb{E}(U) =k\, \beta \qquad\text{et}\qquad Var(U) = k \,\beta^2$$

si $U\sim Gamma(k,\beta)$.

Question 1)¶

A l'aide des formules de l'epérance et de la variance déterminer des estimateurs consistant de $k^\ast$ et $\beta^\ast$.

On pourra comparer son résultat avec celui renvoyé par la commande fit du module scipy.stats qui utilise l'estimateur du maximum de vraisemblance.

In [17]:
sps.gamma.fit(X, floc=0)
Out[17]:
(1.5158713238266452, 0, 1.997800186454287)

Question 2)¶

Déterminer les intervalles de confiance Bootstrap pour $k^\ast$ et $\beta^\ast$ obtenus avec les estimateurs de la question 1) et ceux obtenus avec l'estimateur du maximum de vraisemblance. Les comparer.

Question 3)¶

Comparer les distributions Bootstrap obtenues avec les estimateurs de la question 1 et celles obtenues avec l'estimateur du maximum de vraisemblance. On pourra prendre $B=10000$. Quel estimateur vous paraît le meilleur?

Question 4)¶

Evaluer le biais Bootstrap, la variance Bootstrap, et l'erreur moyenne quadratique Bootstrap des estimateurs de la questions 1 et les comparer avec ceux du maximum de vraisemblance. Que peut-on en conclure?

Solutions:¶

Question 1)¶

In [18]:
k_hat = X.mean()**2 / X.var()
k_hat
Out[18]:
1.5655993534860209
In [19]:
beta_hat = X.var() / X.mean()
beta_hat
Out[19]:
1.934344190062684

Les résultats sont assez similaire à ceux du maximum de vraisemblance sans être exactement les mêmes.

Question 2)¶

Commençons par générer les échantillons Bootstrap.

In [20]:
def dist_boot_gamma(data, B):
    
    k_hat = data.mean()**2 / data.var()
    beta_hat = data.var() / data.mean()

    n = len(data)
    
    X = sps.gamma(k_hat, scale=beta_hat).rvs((B, n))
    
    k_boot = X.mean(axis=1)**2 / X.var(axis=1)
    beta_boot = X.var(axis=1) / X.mean(axis=1)
    
    return k_boot, beta_boot
In [21]:
B = 20000
k_boot, beta_boot = dist_boot_gamma(X, B)
k_tilde, beta_tilde = np.sort(k_boot), np.sort(beta_boot)

On peut aussi obtenir ceux générés avec l'estimateur du maximum de vraisemblance

In [22]:
def dist_boot_gamma_MLE(data, B):
    
    k_hat = data.mean()**2 / data.var()
    beta_hat = data.var() / data.mean()

    n = len(data)
    k_boot, beta_boot = np.zeros(B), np.zeros(B)
    
    X = sps.gamma(k_hat, scale=beta_hat).rvs((B, n))
    
    for j in range(B):
        k_boot[j], _, beta_boot[j] = sps.gamma.fit(X[j,:], floc=0)
    
    return k_boot, beta_boot
In [23]:
k_boot_MLE, beta_boot_MLE = dist_boot_gamma_MLE(X, B)
k_tilde_MLE, beta_tilde_MLE = np.sort(k_boot_MLE), np.sort(beta_boot_MLE)
In [24]:
alpha = 0.05

IC1_k = k_tilde[int(np.floor(alpha * B/2))]
IC2_k = k_tilde[int(np.floor((1-alpha/2) * B))]

IC1_beta = beta_tilde[int(np.floor(alpha * B/2))]
IC2_beta = beta_tilde[int(np.floor((1-alpha/2) * B))]

print("les bornes de l'IC pour k sont %s et %s" % (IC1_k, IC2_k))
print("les bornes de l'IC pour beta sont %s et %s" % (IC1_beta, IC2_beta))

les bornes de l'IC pour k sont 1.4464576277946313 et 1.695529149386269
les bornes de l'IC pour beta sont 1.7700821133642024 et 2.110995617624618
In [25]:
IC1_k_MLE = k_tilde_MLE[int(np.floor(alpha * B/2))]
IC2_k_MLE = k_tilde_MLE[int(np.floor((1-alpha/2) * B))]

IC1_beta_MLE = beta_tilde_MLE[int(np.floor(alpha * B/2))]
IC2_beta_MLE = beta_tilde_MLE[int(np.floor((1-alpha/2) * B))]

print("les bornes de l'IC pour k sont %s et %s" % (IC1_k_MLE, IC2_k_MLE))
print("les bornes de l'IC pour beta sont %s et %s" % (IC1_beta_MLE, IC2_beta_MLE))

les bornes de l'IC pour k sont 1.4805872892189715 et 1.6592967578911497
les bornes de l'IC pour beta sont 1.806081141092021 et 2.0621232034262773

Question 3)¶

In [26]:
figsize(9, 9)

plt.subplot(221)

plt.hist(k_boot, bins=80, density=True, alpha=0.2, label="échantillon");
plt.hist(k_boot_MLE, bins=80, density=True, alpha=0.2, label="échantillon MLE");
plt.vlines(k, 0, 4, label="vraie valeur")
plt.title("histogramme échantillon bootstrap pour $k*$")
plt.legend()

plt.subplot(222)

y = np.linspace(1/B, 1, B)
plt.plot(k_tilde, y)
plt.plot(k_tilde_MLE, y)
plt.title("fonction de répartition bootstrap pour $k*$")

plt.subplot(223)

plt.hist(beta_boot, bins=80, density=True, alpha=0.2, label="échantillon");
plt.hist(beta_boot_MLE, bins=80, density=True, alpha=0.2, label="échantillon MLE");
plt.vlines(beta, 0, 2, label="vraie valeur")
plt.title("histogramme échantillon bootstrap pour beta*")
plt.legend()

plt.subplot(224)

y = np.linspace(1/B, 1, B)
plt.plot(beta_tilde, y)
plt.plot(beta_tilde_MLE, y)
plt.title("fonction de répartition bootstrap pour beta*");
No description has been provided for this image

L'estimateur du maximum de vraisemblance paraît ici meilleur, car les distributions Bootstrap associées semblent légérement plus centrées sur les vraies valeurs, il semble y avoir moins de variance (moins de dispersion).

Question 4)¶

Commençons par l'estimation des biais.

In [27]:
k_hat_MLE, _, beta_hat_MLE = sps.gamma.fit(X, floc=0)
In [28]:
biais_boot = k_boot.mean() - k_hat
biais_boot
Out[28]:
0.004168818711545175
In [29]:
biais_boot_MLE = k_boot_MLE.mean() - k_hat_MLE
biais_boot_MLE
Out[29]:
0.05177565556255437

Pour les variances

In [30]:
var_boot = k_boot.var()
var_boot
Out[30]:
0.004018804633881902
In [31]:
var_boot_MLE = k_boot_MLE.var()
var_boot_MLE
Out[31]:
0.0020597578971435497

On utilise le fait que l'erreur quadratique moyenne (RMSE-Root Mean Square Error) Bootstrap est donnée par

$$ RMSE = \widehat{biais}^2 + \widehat{Var} $$

In [32]:
RMSE = biais_boot**2 + var_boot 
RMSE
Out[32]:
0.004036183683331631
In [33]:
RMSE_MLE = biais_boot_MLE**2 + var_boot_MLE 
RMSE_MLE
Out[33]:
0.004740476406075817

La méthode Bootstrap paramétrique avec la méthode d'estimation par maximum de vraisembalance semble être meilleurs, du point de vue de l'erreur quadratique moyenne, que celle utilisant la méthode des moments pour l'estimation de $k$.

Pour l'estimation de $\beta$, on procède de la même manière.

In [34]:
biais_boot = beta_boot.mean() - beta_hat
biais_boot
Out[34]:
-0.002114251043044524
In [35]:
biais_boot_MLE = beta_boot_MLE.mean() - beta_hat_MLE
biais_boot_MLE
Out[35]:
-0.06469745253641324
In [36]:
var_boot = beta_boot.var()
var_boot
Out[36]:
0.007381955898589861
In [37]:
var_boot_MLE = beta_boot_MLE.var()
var_boot_MLE
Out[37]:
0.004284551352523174
In [38]:
RMSE = biais_boot**2 + var_boot 
RMSE
Out[38]:
0.007386425956062876
In [39]:
var_boot_MLE = k_boot_MLE.var()
var_boot_MLE
Out[39]:
0.0020597578971435497

La méthode Bootstrap paramétrique avec la méthode d'estimation par maximum de vraisembalance semble être aussi meilleurs, du point de vue de l'erreur quadratique moyenne, que celle utilisant la méthode des moments pour l'estimation de $\beta$. Ceci confirme donc l'impression qu'on avait sur les graphiques de la question 3.