Rééchantillonnez avec le bootstrap et les permutations

Votre modèle obtient une AUC de 0,84 sur le jeu de test. Quelle est la précision de ce chiffre ? Avec un autre jeu de test de même taille, auriez-vous obtenu 0,83 ou 0,75 ? Pour une moyenne, une formule donne l’intervalle de confiance. Pour une médiane, une AUC, un F1 ou un écart de taux de conversion, les formules sont compliquées, approximatives ou inexistantes.

Les méthodes de rééchantillonnage contournent le problème en remplaçant les calculs théoriques par de la simulation à partir des données elles-mêmes. Ce chapitre présente les deux principales : le bootstrap, pour estimer l’incertitude d’une statistique, et le test de permutation, pour tester une hypothèse. Ce sont sans doute les outils de ce cours que vous utiliserez le plus souvent en pratique.

Le bootstrap : l’échantillon comme substitut de la population

L’idée

Si l’on pouvait tirer autant d’échantillons que l’on veut dans la population, on calculerait la statistique sur chacun et on verrait directement comment elle varie. Mais on ne dispose que d’un échantillon. L’idée du bootstrap, proposée par Bradley Efron en 1979, est de faire comme si l’échantillon était la population : on remplace la loi inconnue $F$ par la fonction de répartition empirique $F_n$ (chapitre 1), et on tire dans $F_n$.

Tirer un échantillon de taille $n$ dans $F_n$, c’est tirer $n$ fois avec remise parmi les $n$ observations. Certaines observations apparaissent plusieurs fois, d’autres pas du tout.

Chaque échantillon bootstrap mélange doublons et absents. Selon que la valeur 2 400 est tirée ou non, la moyenne passe de 130 à plus de 320 ms, alors que la médiane bouge à peine.

L’algorithme :

1. calculer la statistique $\hat\theta$ sur l’échantillon d’origine ;
2. répéter $B$ fois (typiquement $B = 10\,000$) : tirer un échantillon bootstrap de taille $n$ avec remise, et y calculer la statistique $\hat\theta^*_b$ ;
3. utiliser la distribution des $B$ valeurs $\hat\theta^*_1, \dots, \hat\theta^*_B$ comme approximation de la distribution de $\hat\theta$ : son écart-type estime l’erreur type de $\hat\theta$, ses quantiles donnent un intervalle de confiance.

L’intervalle percentile à 95 % est simplement formé des quantiles 2,5 % et 97,5 % des $\hat\theta^*_b$.

Le fil rouge : la latence médiane

Reprenons les douze temps de chargement du chapitre 1 (en ms) : 120, 135, 128, 142, 131, 125, 138, 129, 2 400, 133, 127, 140. La médiane vaut 132 ms. Quel intervalle de confiance lui donner ?

Exemple

Avec $B = 10\,000$ échantillons bootstrap, l’intervalle percentile à 95 % de la médiane est d’environ [127,5 ; 139] ms. Celui de la moyenne (320,7 ms) est d’environ [129 ; 700] ms.

Pourquoi un intervalle aussi large pour la moyenne ? Un échantillon bootstrap contient au moins une fois la valeur 2 400 avec une probabilité $1 - (11/12)^{12} \approx 0{,}65$, et parfois deux ou trois fois. La moyenne bootstrap saute donc entre environ 130 ms (sans la valeur extrême) et plusieurs centaines de millisecondes. Le bootstrap rend visible l’instabilité de la moyenne face à une valeur extrême, là où la formule classique $\bar x \pm 1{,}96\, s/\sqrt n$ la masquerait derrière une hypothèse de normalité.

Les quantiles 2,5 % et 97,5 % (pointillés) donnent l'intervalle percentile. Celui de la moyenne est énorme parce que sa distribution bootstrap saute d'un paquet à l'autre selon le nombre de fois où 2 400 est tirée.

Les valeurs exactes des bornes varient légèrement d’une exécution à l’autre, puisqu’elles reposent sur des tirages aléatoires ; fixez la graine du générateur pour des résultats reproductibles.

import numpy as np
from scipy import stats

lat = np.array([120, 135, 128, 142, 131, 125, 138, 129, 2400, 133, 127, 140])
rng = np.random.default_rng(42)

# À la main : B échantillons avec remise
B = 10_000
med_boot = np.array([np.median(rng.choice(lat, size=lat.size, replace=True)) for _ in range(B)])
print(np.percentile(med_boot, [2.5, 97.5]))     # environ [127.5 139.]

# Avec scipy (les données sont passées sous forme de tuple)
res = stats.bootstrap((lat,), np.median, n_resamples=B, method="percentile", random_state=rng)
print(res.confidence_interval)

En R, le paquet boot demande une fonction qui reçoit les données et un vecteur d’indices :

library(boot)
lat <- c(120, 135, 128, 142, 131, 125, 138, 129, 2400, 133, 127, 140)
b <- boot(lat, statistic = function(d, i) median(d[i]), R = 10000)
boot.ci(b, type = c("perc", "bca"))

Quel intervalle choisir ?

Intervalle	Principe	Quand l’utiliser
Percentile	quantiles 2,5 % et 97,5 % des $\hat\theta^*$	simple et intuitif, correct si la distribution bootstrap est à peu près symétrique et sans biais
Basique (ou pivot)	$[2\hat\theta - q_{97,5}\,;\ 2\hat\theta - q_{2,5}]$	corrige un décalage, mais peut sortir des valeurs possibles
BCa (bias-corrected and accelerated)	percentile corrigé du biais et de l’asymétrie	le meilleur choix général ; défaut de scipy.stats.bootstrap

À retenir

Le bootstrap estime la variabilité d’une statistique en la recalculant sur des échantillons tirés avec remise dans les données. Il s’applique à presque n’importe quelle statistique (médiane, AUC, F1, rapport de deux moyennes) sans formule théorique.

Les limites du bootstrap

Le bootstrap n’est pas magique. Il repose sur l’idée que $F_n$ ressemble à $F$, ce qui échoue dans plusieurs cas :

• très petits échantillons : avec 5 observations, $F_n$ est une image grossière de $F$, et l’intervalle est trop étroit ;
• statistiques extrêmes : le maximum d’un échantillon bootstrap ne peut jamais dépasser le maximum observé, donc le bootstrap du maximum (ou d’un quantile à 99,9 %) est faux ;
• données dépendantes : séries temporelles, mesures répétées sur un même client. Tirer les observations une à une casse la dépendance ; il faut tirer des blocs (block bootstrap) ou des grappes entières (tous les enregistrements d’un même client ensemble).

Piège

Le rééchantillonnage doit reproduire la façon dont les données ont été produites. Si votre jeu de test contient plusieurs images par patient, rééchantillonnez les patients, pas les images. Sinon, l’intervalle sera bien trop optimiste.

Le test de permutation : fabriquer la loi sous H₀

L’idée

Revenons au test A/B du chapitre 3 : six sessions de la version A, sept de la version B. Sous l’hypothèse nulle « les deux versions produisent la même loi de durées », l’étiquette A ou B d’une session n’a aucune importance : on dit que les étiquettes sont échangeables. On peut donc les redistribuer au hasard, recalculer la statistique, et voir à quel point la valeur observée est exceptionnelle parmi toutes ces redistributions.

L’algorithme :

1. calculer la statistique observée $T_{\text{obs}}$, par exemple l’écart des moyennes $\bar y_B - \bar y_A$ ;
2. pour chaque façon de répartir les $N$ observations en deux groupes de tailles $n_1$ et $n_2$ (ou pour $B$ répartitions tirées au hasard), recalculer $T^*$ ;
3. la p-valeur est la proportion des $T^*$ au moins aussi extrêmes que $T_{\text{obs}}$.

Exemple

Les moyennes valent 4,35 min pour A et 7,81 min pour B : $T_{\text{obs}} \approx 3{,}46$ min. Il y a $\binom{13}{6} = 1\,716$ façons de répartir les 13 sessions en un groupe de 6 et un groupe de 7. Parmi elles, 21 donnent un écart $T^* \ge 3{,}46$ : la p-valeur unilatérale vaut $21/1\,716 \approx 0{,}012$. En bilatéral, 46 répartitions donnent un écart au moins aussi grand en valeur absolue, soit $p \approx 0{,}027$.

C’est très proche du résultat de Mann-Whitney ($p \approx 0{,}022$). Et pour cause : Mann-Whitney est un test de permutation, appliqué aux rangs plutôt qu’aux valeurs. C’est d’ailleurs ainsi que nous avions obtenu sa loi exacte.

La loi de permutation est la loi de T sous H₀, fabriquée à partir des données. La valeur observée tombe dans la queue : 21 répartitions à droite, 25 à gauche.

Quand le nombre de répartitions est trop grand pour les énumérer (plus de 184 000 pour deux groupes de 10, et astronomique au-delà), on en tire $B$ au hasard. On calcule alors la p-valeur par

$$p = \frac{1 + \#\lbrace b : |T^*_b| \ge |T_{\text{obs}}| \rbrace}{B + 1}$$

où l’on compte les permutations au moins aussi extrêmes que l’observation. Le « +1 » au numérateur et au dénominateur compte l’échantillon observé lui-même parmi les permutations : il garantit que $p$ n’est jamais nul et que le test respecte son niveau.

Données appariées : inverser les signes

Pour des données appariées (chapitre 2), on ne permute pas entre groupes : sous $H_0$, chaque différence a autant de chances d’être positive que négative. On inverse les signes au hasard. Pour les dix pages du cache, sur les $2^{10} = 1\,024$ inversions possibles, 12 donnent une somme des différences au moins aussi grande en valeur absolue que la somme observée (−145 ms), d’où $p \approx 0{,}012$.

Sous H₀, chaque signe est à pile ou face. La somme observée, −145 ms, fait partie des 12 cas extrêmes sur 1 024 : p ≈ 0,012.

import numpy as np
from scipy import stats

A = np.array([3.2, 5.1, 2.4, 7.8, 4.0, 3.6])
B = np.array([6.3, 8.9, 4.4, 12.5, 7.1, 5.7, 9.8])

def ecart_moyennes(x, y):
    return np.mean(y) - np.mean(x)

res = stats.permutation_test((A, B), ecart_moyennes,
                             permutation_type="independent",
                             alternative="two-sided", n_resamples=9999)
print(res.statistic, res.pvalue)   # 3.46 ; p ≈ 0.024 (1 716 répartitions < 9 999 : test exact)

Deux remarques sur ce code. D’abord, comme le nombre de répartitions (1 716) est inférieur à n_resamples, scipy les énumère toutes : le test est exact. Ensuite, en bilatéral, scipy double la plus petite des deux p-valeurs unilatérales ($2 \times 21/1\,716 \approx 0{,}024$), au lieu de compter les écarts en valeur absolue (0,027). Les deux conventions sont valables ; elles diffèrent quand la loi de permutation n’est pas symétrique. Pour des données appariées, utilisez permutation_type="samples".

En R, une boucle suffit :

A <- c(3.2, 5.1, 2.4, 7.8, 4.0, 3.6)
B <- c(6.3, 8.9, 4.4, 12.5, 7.1, 5.7, 9.8)
obs <- mean(B) - mean(A)
x <- c(A, B)
set.seed(1)
perm <- replicate(9999, { idx <- sample(13, 7); mean(x[idx]) - mean(x[-idx]) })
(1 + sum(abs(perm) >= abs(obs))) / (9999 + 1)    # 0.0237 avec cette graine (valeur exacte : 46/1716 = 0.027)

Attention

Bootstrap et permutation ne répondent pas à la même question. Le bootstrap tire avec remise et estime la variabilité d’une statistique (intervalle de confiance). La permutation tire sans remise en redistribuant les étiquettes, et simule la loi de la statistique sous $H_0$ (p-valeur). Les confondre donne des p-valeurs fausses.

Application : évaluer des modèles

Un intervalle de confiance pour une AUC

Pour chiffrer l’incertitude d’une AUC mesurée sur un jeu de test, on rééchantillonne les exemples du jeu de test avec remise, on recalcule l’AUC à chaque fois, puis on prend les quantiles. Deux précautions : tirer séparément dans les positifs et dans les négatifs (bootstrap stratifié), pour que chaque échantillon contienne les deux classes, et respecter les grappes éventuelles (plusieurs exemples par patient ou par client).

import numpy as np
from scipy import stats

def auc(scores_pos, scores_neg):
    # AUC = U / (n+ n-), calculée par les rangs (chapitre 3)
    U = stats.mannwhitneyu(scores_pos, scores_neg).statistic
    return U / (len(scores_pos) * len(scores_neg))

res = stats.bootstrap((scores_pos, scores_neg), auc, n_resamples=5000,
                      paired=False, vectorized=False, random_state=0)
print(res.confidence_interval)   # IC à 95 % (BCa par défaut)

Ici, scores_pos et scores_neg sont les scores attribués par le modèle aux exemples positifs et négatifs du jeu de test. En passant les deux échantillons séparément, scipy rééchantillonne chacun dans sa propre classe, ce qui réalise le bootstrap stratifié.

Comparer deux modèles sur le même jeu de test

Deux modèles ont été évalués sur les mêmes exemples de test : leurs résultats sont appariés. Un test de permutation adapté échange, pour chaque exemple et avec une chance sur deux, les prédictions des deux modèles, puis recalcule l’écart de métrique. Sous $H_0$ (les deux modèles se valent), cet échange ne change rien à la loi de l’écart. C’est l’analogue, pour une métrique quelconque, de l’inversion des signes vue plus haut.

En entretien

« Pour un intervalle de confiance sur une AUC ou un F1, je fais un bootstrap sur le jeu de test, stratifié par classe et par client s’il y a des mesures répétées. Pour comparer deux modèles sur le même jeu de test, je fais un test de permutation apparié : j’échange au hasard les prédictions des deux modèles exemple par exemple. »

En résumé

• Le bootstrap remplace la population par l’échantillon : on tire $B$ échantillons avec remise, on recalcule la statistique, et ses quantiles donnent un intervalle de confiance (percentile, basique ou BCa).
• Il s’applique à presque toute statistique, mais échoue pour les très petits échantillons, les statistiques extrêmes et les données dépendantes (il faut alors rééchantillonner par blocs ou par grappes).
• Le test de permutation redistribue les étiquettes, sans remise, pour simuler la loi de la statistique sous $H_0$ ; il est exact quand on énumère toutes les permutations.
• En Monte-Carlo, $p = (1 + \#\lbrace |T^*| \ge |T_{\text{obs}}| \rbrace)/(B + 1)$ ; pour des données appariées, on inverse les signes.
• Mann-Whitney et Wilcoxon sont des tests de permutation appliqués aux rangs.
• En apprentissage automatique : bootstrap stratifié pour l’IC d’une métrique, permutation appariée pour comparer deux modèles sur le même jeu de test.

Exercices

Exercice · Compter les permutations

facile

1. Combien y a-t-il de façons de répartir 10 observations en deux groupes de 5 ? Et 20 observations en deux groupes de 10 ?
2. Pour un test apparié sur 15 paires, combien d’inversions de signes faut-il énumérer ?
3. Dans quels cas faut-il passer à un tirage aléatoire de $B$ permutations ?

Voir la solution

1. $\binom{10}{5} = 252$ et $\binom{20}{10} = 184\,756$.
2. $2^{15} = 32\,768$.
3. Les trois cas restent énumérables par ordinateur. Mais le nombre explose vite : $\binom{40}{20} \approx 1{,}4 \times 10^{11}$, $2^{50} \approx 10^{15}$. Au-delà de quelques centaines de milliers de cas, on tire $B$ permutations au hasard (typiquement 10 000) et on utilise $p = (1 + \#)/(B + 1)$. L’erreur de Monte-Carlo sur $p$ est d’environ $\sqrt{p(1-p)/B}$, soit 0,002 pour $p = 0{,}05$ et $B = 10\,000$.

En code

Python

from math import comb, sqrt

print(comb(10, 5), comb(20, 10))     # 252 184756
print(2 ** 15)                       # 32768
print(comb(40, 20), 2 ** 50)         # 137846528820 (1.4e11) 1125899906842624 (1.1e15)

# Erreur de Monte-Carlo sur p, pour p = 0.05 et B = 10 000 tirages
p, B = 0.05, 10_000
print(sqrt(p * (1 - p) / B))         # 0.0022

R

choose(10, 5); choose(20, 10)        # 252 et 184756
2^15                                 # 32768
choose(40, 20); 2^50                 # 137846528820 (1.4e11) et 1.126e+15

# Erreur de Monte-Carlo sur p, pour p = 0.05 et B = 10 000 tirages
p <- 0.05; B <- 10000
sqrt(p * (1 - p) / B)                # 0.0022

Exercice · Les observations absentes d'un échantillon bootstrap

moyen

1. Pour un échantillon de taille $n$, quelle est la probabilité qu’une observation donnée n’apparaisse pas dans un échantillon bootstrap ?
2. Que vaut cette probabilité quand $n$ devient grand ?
3. Quel lien avec les forêts aléatoires ?

Voir la solution

1. À chacun des $n$ tirages, l’observation a une probabilité $1 - 1/n$ de ne pas être choisie. Les tirages étant indépendants, la probabilité qu’elle ne soit jamais choisie est $(1 - 1/n)^n$.
2. Elle tend vers $e^{-1} \approx 0{,}368$. Un échantillon bootstrap contient donc en moyenne environ 63,2 % des observations distinctes ; les autres sont des doublons. Pour $n = 12$, on trouve $(11/12)^{12} \approx 0{,}352$, ce qui explique que la valeur 2 400 de notre exemple figure dans environ 65 % des échantillons bootstrap.
3. Chaque arbre d’une forêt aléatoire est entraîné sur un échantillon bootstrap (c’est le bagging). Les quelque 37 % d’observations absentes, dites hors sac (out-of-bag), servent de jeu de validation gratuit pour cet arbre : c’est le score OOB.

En code

Python

import numpy as np

# Probabilité qu'une observation donnée soit absente d'un échantillon bootstrap
for n in [12, 100, 10_000]:
    print(n, (1 - 1 / n) ** n)       # 0.352, 0.366, 0.368
print(np.exp(-1))                    # limite : 0.368

# Vérification par simulation : proportion d'observations distinctes tirées
rng = np.random.default_rng(0)
n = 1000
distinctes = [len(np.unique(rng.integers(0, n, size=n))) / n for _ in range(2000)]
print(np.mean(distinctes))           # environ 0.632

R

# Probabilité qu'une observation donnée soit absente d'un échantillon bootstrap
n <- c(12, 100, 10000)
(1 - 1 / n)^n                        # 0.352 0.366 0.368
exp(-1)                              # limite : 0.368

# Vérification par simulation : proportion d'observations distinctes tirées
set.seed(0)
distinctes <- replicate(2000, length(unique(sample(1000, replace = TRUE))) / 1000)
mean(distinctes)                     # environ 0.632

Exercice · Quand le bootstrap échoue

difficile

On veut un intervalle de confiance pour le maximum $\theta$ d’une loi uniforme sur $[0, \theta]$, à partir de $n = 50$ observations, en prenant le maximum de l’échantillon comme estimateur.

1. Quelle est la probabilité qu’un échantillon bootstrap ait exactement le même maximum que l’échantillon d’origine ?
2. Pourquoi l’intervalle percentile est-il alors inutilisable ?
3. Que proposeriez-vous ?

Voir la solution

1. Le maximum bootstrap est égal au maximum observé dès que la plus grande observation est tirée au moins une fois : probabilité $1 - (1 - 1/50)^{50} \approx 1 - 0{,}364 = 0{,}636$.
2. La distribution bootstrap du maximum a donc une masse de près de 64 % sur une seule valeur, et elle ne peut jamais dépasser le maximum observé. Or le vrai $\theta$ est presque toujours au-dessus du maximum observé. L’intervalle percentile, entièrement situé sous le maximum observé, ne contient pratiquement jamais $\theta$ : sa couverture réelle est très loin de 95 %.
3. Utiliser la théorie : pour une loi uniforme, $P(\max \le t) = (t/\theta)^n$, ce qui donne un intervalle exact $[\max\,;\ \max / 0{,}05^{1/n}]$ à 95 %. Plus généralement, pour les statistiques extrêmes, on se tourne vers un bootstrap paramétrique, un bootstrap « m sur n » (échantillons de taille $m$ plus petite que $n$) ou la théorie des valeurs extrêmes.

En code

Python

import numpy as np

rng = np.random.default_rng(0)
theta, n = 1.0, 50
x = rng.uniform(0, theta, size=n)
m = x.max()

# Bootstrap du maximum
max_boot = np.array([rng.choice(x, size=n).max() for _ in range(10_000)])
print(1 - (1 - 1 / n) ** n, np.mean(max_boot == m))   # 0.636 et 0.633
print(m, np.percentile(max_boot, [2.5, 97.5]))        # 0.997 [0.934 0.997] : l'IC reste sous m < theta

# Intervalle exact [max ; max / 0.05^(1/n)] : il contient bien theta = 1
print(m, m / 0.05 ** (1 / n))                         # 0.997 1.059

R

set.seed(0)
theta <- 1; n <- 50
x <- runif(n, 0, theta)
m <- max(x)

# Bootstrap du maximum
max_boot <- replicate(10000, max(sample(x, replace = TRUE)))
c(1 - (1 - 1 / n)^n, mean(max_boot == m))   # 0.636 et 0.650
c(m, quantile(max_boot, c(0.025, 0.975)))   # 0.992 [0.908 0.992] : l'IC reste sous m < theta

# Intervalle exact [max ; max / 0.05^(1/n)] : il contient bien theta = 1
c(m, m / 0.05^(1 / n))                      # 0.992 1.053