Introduction à bootstrap

Principe

• On s'intéresse au paramètre d'une loi de probabilité $F_\theta$, $f_\theta$, $f(x; \theta)$ ou $\mathbb{P}_\theta$

• On ne possède qu'un échantillon $X_1, \ldots, X_n$ issu de $F_\theta$

• $\widehat{\theta} = t\big(X_1, \ldots, X_n\big)$ est un estimateur de $\theta$

Plusieurs questions liées à la distribution de $\widehat\theta$

• Biais de l'estimateur $\widehat\theta$
• Variance de $\widehat\theta$ (sa précision)
• Intervalle de confiance pour $\theta$

Simuler suivant $F_\theta$ ou $f_\theta$

Exemple : supposons que $X \sim \mathcal N(\theta, 1)$ où $\theta = 2$

# simuler 100 réalisations d'une loi normale  suivant une loi N(0, 1)
#import numpy as np
import numpy as np
from scipy.stats import norm, uniform, expon
import matplotlib.pyplot as plt
import pandas as pd
import statsmodels.api as sm
from sklearn import metrics
import warnings
warnings.filterwarnings("ignore")

#norm_rv.random_state = np.random.RandomState(seed=11)
np.random.seed(11)
x = norm.rvs(loc =2, scale =1, size=100)


# simuler N = 10000 jeu de données
n =  100
N =  10000 
X = np.reshape(norm.rvs(loc =2, scale =1, size=n*N), (N,n))
t = X.mean(axis=1)
t.shape
print(t.mean())
print(t.std())
print(t.var())


# comparer des distribution empirique versus empirique 

xx = np.linspace(norm.ppf(0.01, loc=2, scale = 1/np.sqrt(n)), 
                 norm.ppf(0.99, loc=2, scale = 1/np.sqrt(n)), 100)
fig, ax = plt.subplots(1, 1)
ax.plot(xx, norm.pdf(xx, loc=2, scale = 1/np.sqrt(n)), 'r-', lw=5, alpha=0.6, label='densite')
ax.hist(t, density=True, histtype='stepfilled', alpha=0.2, bins =30)

Une petite parenthèse : générateur de nombres aléatoires

La loi de base que l'on simule est la loi uniforme sur $[0, 1]$, un tel générateur est inclus dans tous les langages de programmation et dans tous les logiciels.

• La méthode de congruence linéaire (D.H Lehmer-Lucas 1948), \begin{equation*} g_{n+1} = a g_n + b \mod{m} \end{equation*} où $a,b,m$ et $g_0$ sont à choisir intialement.

• Ce générateur fabrique une suite d'entiers aléatoires allant de $0$ à $m-1$.

• La suite $\frac{g_n}{m-1}$ fournit des nombres sur l'intervalle $[0, 1]$.

• La suite ainsi produite possède une période inférieure ou égale à $m$.
• Un générateur de ce type, nommé RANDU, fut exploité pendant des années dans toutes sortes de logiciels.
• Les paramètres de RANDU sont $m = 2^{35}$, $a = 3125$, $b = 1$ et $g_0=71$ et sa période est de $2^{35}$

Inversion de la fonction de répartition

On veut générer des nombres aléatoires suivant une loi données par sa fonction de répartition $F$. Comment procéder ?

• Générer $U_1, U_2, \ldots, U_n$ suivant une loi uniforme $U(0,1)$.

• On forme \begin{equation*} \big(X_1, X_2, \ldots, X_n\big) = \big(F^{-1}(U_1), F^{-1}(U_2), \ldots, F^{-1}(U_n)\big) \end{equation*}

On va s'intéresser à la loi exponentielle $\mathcal E\big(\lambda\big)$ où $\lambda = 2$

• Rappelons que $F(x) = 1 - e^{-\lambda x }$

• Simuler $100$ réalisations de la loi $\mathcal E(2)$

• Tracer sur la même figure l'histogramme des réalisations ainsi que la vraie densité de probabilité de la loi $\mathcal E(2)$

x_u =  uniform.rvs(size=1000) 
x_e =  expon.ppf(x_u, scale=1/2)

xx = np.linspace(expon.ppf(0.01, scale = 1/2), 
                 expon.ppf(0.99, scale = 1/2), 100)
fig, ax = plt.subplots(1, 1)
ax.hist(x_e, density=True, histtype='stepfilled', alpha=0.2, bins=50)
ax.plot(xx, expon.pdf(xx, scale = 1/2), 'r-', lw=5, alpha=0.6, label='densite')

Simuler une loi normale : Box-Muller

Soit, $U$ et $V$ deux variables aléatoires uniformes sur $[0, 1]$ indépendantes, les variables \begin{equation*} X = \cos(2\pi U) \sqrt{\big(-2 \log(V)\big)} \quad \text{et} \quad Y = \sin(2\pi U) \sqrt{\big(-2 \log(V)\big)} \end{equation*} sont indépendantes et suivent une loi normale $\mathcal N(0, 1)$.

Alternative : rééchantillonnage

$F_\theta$ est en général inconnue !

Le bootstrap consiste donc à faire une simulation à partir de la loi empirique $F_n$ observée (i.e. l'échantillon au lieu de la vraie loi $F_\theta$, qui est inconnue.

• On tire $B$ échantillons $x_1^*, x_2^*, \ldots, x_B^*$ où chaque $x^*_b$ est constitué en tirant avec remise $n$ valeurs de l'échantillon $X_1, \ldots, X_n$.
• On appelle chaque tirage d'un $x_b^*$ un bootstrap

coronary = pd.read_csv("data/coronary.csv", delimiter=';')
coronary = coronary.drop(coronary.columns[0], axis=1)
#print(coronary.head())

y = np.array(coronary['coron'])
X = np.array(coronary.drop(['coron'], axis=1))
print(X[:4,])
X1 = sm.add_constant(X) # adding a constant
print(X1[:4,])
model = sm.Logit(y, X1)
result = model.fit()
result.summary()
#print(result.params)

def boot_logistic_reg(X, y, size=1):
    """bootstrap pour les coefficients de régression logistique."""

    # Set up array of indices to sample from: inds
    inds = np.arange(X.shape[0])

    # Initialize replicates
    bs_params_reps = np.empty(shape=(size, X.shape[1]))

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds)) # sampling the indices (1d array requirement)
        bs_X, bs_y = X[bs_inds], y[bs_inds]
        bs_model = sm.Logit(bs_y, bs_X)
        bs_result = bs_model.fit(disp=False)
        bs_params_reps[i:] = bs_result.params

    return bs_params_reps

toto = boot_logistic_reg(X1, y, 100)

toto.std(axis=0)

 inds = np.arange(X.shape[0])
 print(inds)   
 bs_inds = np.random.choice(inds, size=len(inds))
 print(bs_inds)   

Courbe ROC

Le code suivant permet de tracer la courbe ROC associée à la régression logistique

probs = model.predict(result.params, X1)
print(probs[:10])
fpr, tpr, thresholds = metrics.roc_curve(y, probs)
roc_auc = metrics.auc(fpr, tpr)

## point le plus pro0h1 du coin (0,1)
d = np.sqrt((fpr - 0)**2 + (tpr - 1)**2)
print('Le seuil optimal est %0.3f' % thresholds[np.argmin(d)])

## on peut maintenant tracer la courbe
plt.plot(fpr, tpr, color='darkorange',
         lw=2, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.plot([fpr[np.argmin(d)]], [tpr[np.argmin(d)]], marker='o', markersize=5, color="red")
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()

1. Écrire une fonction boot_threshold_auc(X, y, size) qui permet de produire size réplications bootstrap du seuil optimal et size réplications de l'aire sous la courbe ROC.

2. Proposer un intervalle de confiance empirique de l'aire sous la courbe ROC