TP1 Méthode de Bootstrap (suite)¶

In [1]:
import numpy as np
import scipy.stats as sps 
import matplotlib.pyplot as plt
from IPython.core.pylabtools import figsize

5) Estimateur du minimax vs estimateur des moments¶

Dans cette partie, on va considérer un échantillon $(X_1,\dots,X_n)$ un échantillon de loi de Bernoulli de parmaètre $p\in]0,1[$. On rappelle que $\mathbb{E}(X_1)=p$ et on va comparer l'estimateur de la moyenne $\bar X_n$ (qui est aussi l'estimateur du maximum de vraisemblance dans ce cas), avec l'estimateur du minimax avec coût quadratique

$$ \hat \theta_n = \frac{\sqrt{n}}{\sqrt{n}+1} \bar X_n +\frac{1}{2(\sqrt{n}+1)}.$$

Question:¶

Comparer les RMSE Bootstrap de ces deux estimateurs en fonction de $p$ et de la taille de l'échantillon $n$. On pourra faire des graphiques avec en abscisse $p$, et on pourra prendre $B=2000$, $n$ dans $\{30,2000\}$ et $p$ entre $0$ et $1$ avec 2000 points (np.linspace(0, 1, 30)). Que peut-on conclure?

In [2]:
def Moment_p_boot(X, p, B):
    
    n = len(X)
    p_hat = X.mean()
    
    p_boot = np.zeros(B)
    
    for b in range(B):
        X = sps.bernoulli(p_hat).rvs(n)
        p_boot[b] = X.mean()
        
    biais_boot = p_boot.mean() - p_hat
    var_boot = p_boot.var()
 
    return biais_boot**2 + var_boot


def MM_p_boot(X, p, B):
    
    n = len(X)
    sqrt_n = np.sqrt(n)
    c1 = sqrt_n / (sqrt_n + 1)
    c2 = 0.5 / (sqrt_n + 1)
    
    p_hat = c1 * X.mean() + c2
    
    p_boot = np.zeros(B)
    
    for b in range(B):
        X = sps.bernoulli(p_hat).rvs(n)
        p_boot[b] = c1 * X.mean() + c2
        
    biais_boot = p_boot.mean() - p_hat
    var_boot = p_boot.var()
    
    return biais_boot**2 + var_boot
In [3]:
def plot_RMSE(P, N, B):
    
    moment_boot = np.zeros((len(N), len(P)))
    mm_boot = np.zeros((len(N), len(P)))
    
    for (i,n) in enumerate(N):
        for (j,p) in enumerate(P):
            X = sps.bernoulli(p).rvs(n)
            moment_boot[i,j] = Moment_p_boot(X, p, B)
            mm_boot[i,j] = MM_p_boot(X, p, B)
            
    return moment_boot, mm_boot
In [4]:
P = np.linspace(0, 1, 20)
N = [30, 2000]
B = 2000

moment_boot, mm_boot = plot_RMSE(P, N, B)
In [5]:
plt.subplot(121)
plt.plot(P, moment_boot[0,:], label="Moment/MLE")
plt.plot(P, mm_boot[0,:], label="MiniMax")
plt.xlabel("p")
plt.ylabel("RMSE")
plt.legend()

plt.subplot(122)
plt.plot(P, moment_boot[1,:], label="Moment/MLE")
plt.plot(P, mm_boot[1,:], label="MiniMax")
plt.xlabel("p")
plt.ylabel("RMSE")
plt.legend()

figsize(15,5)
No description has been provided for this image

En conclusion, l'estimateur du minimax à une RMSE constante en fonction de $p$, mais pas l'estimateur du maximum de vraisemblance. Pour des grandes valeurs de $n$ l'estimateur du maximum de vraisemblance semble être meilleur que celui du minimax. Par contre, pour des petites valeurs de $n$ l'estimateur du minimax semble être plus performant pour des valeurs de $p$ (inconnu en pratique) entre $0.4$ et $0.6$.

6) Bootstrap non paramétrique¶

La méthode de Bootstrap est plus connue sous sa forme non-paramétrique, mais le principe est similaire. Au lieu de resimuler un échantillon avec une certaine loi paramétrique, avec des paramètres estimés, on réutilise l'échantillon original en le réenchantillonant avec remplacement. Autrement dit, les nouveaux échantillons sont crées en tirant des valeurs de manière uniforme dans l'échantillon original avec remplacement.

Pour cela on utilisera la fonction np.random.choice(S, n) où $n$ est la forme de l'échantillon souhaité et $S$ est un ensemble de valeur (dans notre cas l'ensemble des valeurs de notre échantillon). Un exemple:

In [6]:
S = np.arange(1,5)
n = (5,6)
np.random.choice(S, n)
Out[6]:
array([[3, 2, 4, 2, 3, 1],
       [1, 1, 2, 3, 2, 4],
       [3, 3, 3, 3, 3, 2],
       [2, 2, 3, 4, 2, 2],
       [1, 3, 1, 2, 3, 2]])

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

  • créer de nouveaux échantillons $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(B)}$ (bootstrap sample) à partir de l'échantillon original. Pour chaque échantillon, on tire des valeurs de manière uniforme (avec remplacement) dans l'échantillon original $\mathbf{X} = (X_1,\dots,X_n)$.
  • puis à partir de ces échantillons procéder à de nouvelles estimations $\hat \theta^{(1)},\dots,\hat \theta^{(B)}$ (bootstrap replicates)

Reprenons notre exemple d'un échantillon de loi normale. Faisons comme si nous ne savions pas que cela venaint d'une loi normale, et calculons les distributions Bootstrap pour la moyenne et la variance. On comparera cette distribution avec celle obtenue avec l'approche paramétrique.

In [7]:
n = 30
mu = 0.
sig = np.sqrt(2)

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

Les estimations pour cet echantillon sont

In [8]:
mu_hat = X.mean()
sig_hat = X.std()

mu_hat, sig_hat
Out[8]:
(-0.027315054781262706, 1.4426289708897668)

Définissons la fonction qui va nous produire les échantillon Bootstrap nonparamétrique

In [9]:
def dist_boot_nonpara(data, B):
        
    n = len(data)
    
    X = np.random.choice(data,(B, n))
    
    mu_boot = X.mean(axis=1)
    sig_boot = X.std(axis=1)
    
    return mu_boot, sig_boot

et redennons la fonction qui va nous produire les échantillon Bootstrap paramétrique du TP précédent.

In [12]:
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 [14]:
B = 20000

mu_boot, sig_boot = dist_boot(X, B)
mu_boot_nonpara, sig_boot_nonpara = dist_boot_nonpara(X, B)

Jettons un coup d'oeil aux distributions Bootstrap de la moyenne et de la variance

In [15]:
plt.subplot(121)
plt.hist(mu_boot, bins=80, density=True, alpha=0.2, label="param")
plt.hist(mu_boot_nonpara, bins=50, density=True, alpha=0.2, label="non param")
plt.vlines(mu, 0, 2, label="vraie valeur")
plt.legend()

plt.subplot(122)
plt.hist(sig_boot, bins=80, density=True, alpha=0.2, label="param")
plt.hist(sig_boot_nonpara, bins=50, density=True, alpha=0.2, label="non param")
plt.vlines(sig, 0, 2, label="vraie valeur")
plt.legend()

figsize(15,5)
No description has been provided for this image

Alors que pour l'estimation de la moyenne les méthodes paramétrique et non paramétrique semblent être aussi précises, la méthode paramétrique semble plus précise que la non-paramétrique pour l'estimation de la variance.

Calculons les RMSE de la moyenne et de la variance

In [16]:
biais_mu = mu_boot.mean() - mu_hat
var_mu = mu_boot.var()
RMSE_mu = biais_mu**2 + var_mu

biais_mu_nonpara = mu_boot_nonpara.mean() - mu_hat
var_mu_nonpara = mu_boot_nonpara.var()
RMSE_mu_nonpara = biais_mu_nonpara**2 + var_mu_nonpara

RMSE_mu, RMSE_mu_nonpara
Out[16]:
(0.06883268406437261, 0.06957082195076598)

Pour la moyenne les RMSE sont très proches.

In [17]:
biais_sig = sig_boot.mean() - sig_hat
var_sig = sig_boot.var()
RMSE_sig = biais_sig**2 + var_sig

biais_sig_nonpara = sig_boot_nonpara.mean() - sig_hat
var_sig_nonpara = sig_boot_nonpara.var()
RMSE_sig_nonpara = biais_sig_nonpara**2 + var_sig_nonpara

RMSE_sig, RMSE_sig_nonpara
Out[17]:
(0.03625884963949387, 0.03248256382346437)

On voit bien une RMSE légérement plus petite pour la méthode paramétrique ce qui confirme ce que suggère le graphique ci-dessus.

7) Application à l'évolution des températures en France¶

Commençons par charger les données du fichier data_France.csv.

In [18]:
import pandas as pd
In [19]:
data = pd.read_csv("data_France.csv")

Pour obtenir un aperçu des données

In [20]:
data.head()
Out[20]:
Year January February March April May June July August September October November December Winter Spring Summer Fall MeteorologicalYear
0 1961 0.098 3.417 1.726 2.680 -0.249 0.543 -0.511 -0.072 2.906 0.820 -0.313 0.435 0.828 1.386 -0.013 1.138 0.835
1 1962 1.602 -0.588 -2.716 -0.028 -1.209 -0.375 -0.522 0.804 -0.026 0.549 -1.503 -2.711 0.483 -1.318 -0.031 -0.327 -0.298
2 1963 -4.597 -4.475 -0.135 0.425 -0.684 -0.186 0.389 -1.456 -0.494 0.029 2.512 -2.600 -3.928 -0.131 -0.418 0.682 -0.949
3 1964 -1.664 0.972 -0.942 0.477 1.744 0.834 1.324 0.323 1.244 -1.790 0.424 -1.016 -1.097 0.426 0.827 -0.041 0.029
4 1965 0.346 -2.743 0.008 -0.472 -0.097 0.118 -1.277 -0.861 -1.988 0.804 -0.334 1.471 -1.138 -0.187 -0.673 -0.506 -0.626

Pour obtenir le nom des colonnes

In [21]:
data.columns
Out[21]:
Index(['Year', 'January', 'February', 'March', 'April', 'May', 'June', 'July',
       'August', 'September', 'October', 'November', 'December', 'Winter',
       'Spring', 'Summer', 'Fall', 'MeteorologicalYear'],
      dtype='object')

On va s'intéresser à ce qu'il se passe à l'année

In [22]:
df = data[ ["Year", "MeteorologicalYear"] ]
df.head()
Out[22]:
Year MeteorologicalYear
0 1961 0.835
1 1962 -0.298
2 1963 -0.949
3 1964 0.029
4 1965 -0.626

et pour dessiner nos données

In [23]:
df.plot(x = "Year", y="MeteorologicalYear", kind="scatter", alpha=0.5);
No description has been provided for this image
In [24]:
df.plot(x = "Year", y="MeteorologicalYear", 
        marker="o", ylabel = "variation température (C)", alpha=0.5);
No description has been provided for this image

Remarque:¶

La colonne Year du fichier de données prend des valeurs relativement grandes par rapport aux variations de température. Il est en général préférable de "standardiser" ou "normaliser" cette colonne. Nous ne le ferons pas ici par souci de simplification, mais sachez que c'est une bonne pratique de le faire.

On peut remarquer que l'augmentation est similaire à une droite. On va donc modéliser ces données de la façon suivante

$$ Y_j = \alpha + \beta \, T_j + \varepsilon_j,\qquad j \in\{1,\dots,n\}.$$

où

  • $n$ est le nombre d'observation;
  • $Y_j$ les variations de température observées par rapport à la valeur référence (une moyenne des température entre 1960 et 1980);
  • $T_j$ les année d'observations
  • $\varepsilon_j$ des variables aléatoires iid de loi $\mathcal{N}(0,\sigma^2)$;
  • Les paramètres du modèle $\alpha$, $\beta$ et $\sigma$ sont à estimer à partir des données.

Question 1)¶

En remarquant que

$$ Y_j \sim \mathcal{N}(\alpha + \beta \, T_j,\sigma^2)\qquad j \in\{1,\dots,n\} $$

Déterminer la fonction de vraisemblance $\mathcal{L}(Y_1,\dots,Y_n; \alpha,\beta,\sigma)$ associée aux observations $(Y_1,\dots,Y_n)$.

Question 2)¶

Résoudre l'équation

$$\nabla \log \mathcal{L}(Y_1,\dots,Y_n; \alpha,\beta,\sigma) = 0,$$

où le gradiant porte sur les variables $\alpha$, $\beta$ et $\sigma^2$. On doit trouver comme solution

$$ \hat \beta = \frac{cov(Y,T)}{var(T)},\qquad \hat \alpha = \bar Y_n - \hat \beta \, \bar T_n \qquad\text{et}\qquad \hat \sigma^2 = \frac{1}{n}\sum_{j=1}^n (Y_j - \hat\alpha - \hat \beta\,T_j)^2,$$

où

$$ cov(Y,T) = \frac{1}{n}\sum_{j=1}^n (Y_j - \bar Y_n)(T_j - \bar T_n) \qquad\text{et}\qquad var(T) = \frac{1}{n}\sum_{j=1}^n (Y_j - \bar Y_n)^2.$$

Ceci nous donne les estimateurs du maximum de vraisemblance, autrement dit la solution de $$ \max_{\alpha,\beta,\sigma} \log \mathcal{L}(Y_1,\dots,Y_n; \alpha,\beta,\sigma).$$

Question 3)¶

Calculer une estimation de $\hat \alpha$, $\hat \beta$, $\hat \sigma^2$ et tracer la droite d'équation $t\mapsto \hat \alpha + \hat \beta \,t$ superposée aux données. On pourra utiliser la fonction cov de numpy qui renvoit la matrice de covariance (empirique) associée à deux jeux de données.

Question 4)¶

Déterminer la distribution Bootstrap des paramètres $\alpha$, $\beta$, et $\sigma$. On pourra prendre $B=2000$, et illustrer leurs histogrammes normalisés.

Question 5)¶

Tracer les droites correspondantes aux valeurs de l'échantillon Bootstrap et superposer y les données.

Question 6)¶

Comment pourrait-on évaluer la probabilité qu'il y ait une augmentation des températures moyenne sur l'année à partir de l'échantillon Bootstrap?

Question 7)¶

Avec ce modèle de droite, donner une estimation de l'année pour laquelle on dépasse les $3^o C$ d'augmantation des températures avec $5\%$ de chance de se tromper.

Solutions¶

question 1)¶

La fonction de vraisemblance est donnée par

$$ \mathcal{L}(Y_1,\dots,Y_n; \alpha,\beta,\sigma) = \prod_{j=1}^n \frac{1}{\sqrt{2\pi \sigma^2}} e^{-(Y_j-\alpha-\beta T_j)^2/(2\sigma^2)} = \frac{1}{(2\pi)^{n/2} \sigma^n} e^{-\frac{1}{2\sigma^2}\sum_{j=1}^n (Y_j-\alpha-\beta T_j)^2}$$

question 2)¶

C'est du calcul.

question 3)¶

On calcule ainsi les estimations

In [25]:
cov_mat = np.cov(df.MeteorologicalYear, df.Year)
beta_hat = cov_mat[0,1]/cov_mat[1,1]
alpha_hat = df.MeteorologicalYear.mean() - beta_hat * df.Year.mean()
sigma_hat = np.std( df.MeteorologicalYear - alpha_hat - beta_hat * df.Year )
In [26]:
alpha_hat, beta_hat, sigma_hat
Out[26]:
(-68.70923670368205, 0.0348572180011689, 0.4515710612882542)
In [30]:
df.plot(x = "Year", y="MeteorologicalYear", 
        marker="o", ylabel = "variation température (C)", alpha=0.5, label="data")

years = df.Year.to_numpy()
plt.plot(years, alpha_hat + beta_hat * years, label="linear_approx")
plt.legend();
No description has been provided for this image

question 4)¶

In [31]:
def Estim_boot(df, B):
    
    n = df.shape[0]
    T_mean = df.Year.mean()
    
    cov_mat = np.cov(df.MeteorologicalYear, df.Year)
    beta_hat = cov_mat[0,1]/cov_mat[1,1]
    alpha_hat = df.MeteorologicalYear.mean() - beta_hat * T_mean
    sigma_hat = np.std( df.MeteorologicalYear - alpha_hat - beta_hat * df.Year )
    
    alpha_boot = np.zeros(B)
    beta_boot = np.zeros(B)
    sigma_boot = np.zeros(B)
    
    
    for b in range(B):
        
        Y = alpha_hat + beta_hat * df.Year + sigma_hat * sps.norm().rvs(n)
        
        cov_mat = np.cov(Y, df.Year)
        
        beta_tmp = cov_mat[0,1]/cov_mat[1,1]
        alpha_tmp = Y.mean() - beta_tmp * T_mean
        sigma_tmp = np.std( Y - alpha_tmp - beta_tmp * df.Year )
        
        beta_boot[b] = beta_tmp
        alpha_boot[b] = alpha_tmp
        sigma_boot[b] = sigma_tmp
 
    return alpha_boot, beta_boot, sigma_boot
In [32]:
B= 2000
alpha_boot, beta_boot, sigma_boot = Estim_boot(df, B)
In [33]:
plt.hist(alpha_boot, bins=50, alpha=0.2, density=True);
No description has been provided for this image
In [34]:
plt.hist(beta_boot, bins=50, alpha=0.2, density=True);
No description has been provided for this image
In [35]:
plt.hist(sigma_boot, bins=50, alpha=0.2, density=True);
No description has been provided for this image

question 5)¶

Cette figure nous permet d'évaluer la qualité du modèle (ici une droite) par rapport aux données.

In [37]:
df.plot(x = "Year", y="MeteorologicalYear", 
        marker="o", ylabel = "variation température (C)", alpha=0.5, label="data")

years = df.Year.to_numpy()
for b in range(B):
    plt.plot(years, alpha_boot[b] + beta_boot[b] * years, color="orange", alpha=0.05)
No description has been provided for this image

question 6)¶

On cherche à évaluer $\mathbb{P}(\beta > 0)$ qu'on approche dans notre contexte par $\mathbb{P}(\hat \beta > 0)$ et donc par

$$ \frac{1}{B} \sum_{j=1}^B \mathbf{1}_{\big(\hat \beta^{(j)}>0\big)} \simeq \mathbb{P}(\hat \beta > 0).$$

On a alors

In [38]:
np.mean(beta_boot > 0)
Out[38]:
1.0

question 7)¶

Ici on cherche $T_0$ pour lequel $\alpha + \beta T_0 +\sigma * \varepsilon = 3$, c'est à dire

$$ T_0 = \frac{3-\alpha - \sigma * \varepsilon}{\beta}.$$

où $\varepsilon\sim \mathcal{N}(0,1)$. On peut alors obtenir une distribution Bootstrap pour $T_0$, et obtenir un intervalle de confiance.

Remarque: Comme on est en 2021, il faudra conditionner au fait que $T_0>2021$

In [39]:
T0_boot = (3-alpha_boot - sigma_boot * sps.norm().rvs(B))/beta_boot
T0_boot = T0_boot[T0_boot>2021]

En faisant se filtrage, on obtient bien la loi Bootstrap conditionnellement à $T_0>2021$ mais on perd quelque valeurs dans l'échantillon et donc en précision pour l'estimation.

In [40]:
plt.hist(T0_boot, bins=50, alpha=0.2);
No description has been provided for this image
In [41]:
alpha = 0.025
T0_tilde = np.sort(T0_boot)

IC1_T0 = T0_tilde[int(np.floor(alpha*B/2))]
IC2_T0 = T0_tilde[int(np.floor((1-alpha/2)*B))]

print("les bornes de l'IC pour mu sont %s et %s" % (IC1_T0, IC2_T0))

les bornes de l'IC pour mu sont 2029.8766316166264 et 2095.119509581521