Estimez et comparez par rééchantillonnage

Vous savez maintenant que le vrai risque $R(\hat{f})$ est ce qui compte, que le risque empirique d’apprentissage le sous-estime, et que la capacité du modèle doit être réglée ni trop bas ni trop haut. Reste un problème très concret : $P_{XY}$ est inconnue, donc $R(\hat{f})$ ne se calcule pas. Il faut une stratégie pour l’approcher.

Cette stratégie, c’est le rééchantillonnage (resampling) : on réutilise intelligemment les données dont on dispose pour simuler des « données nouvelles ». Ce chapitre couvre les étapes 6 à 8 de la roue : vérifier la machine, comparer au moins deux approches et dire avec quelle confiance on annonce un résultat. Fil rouge : notre service carrières hésite entre une droite et une méthode non paramétrique pour prédire les salaires, et veut savoir laquelle annoncer aux étudiants.

Réserver des données : l’estimateur du test

L’idée la plus simple consiste à couper l’échantillon en deux : une part pour apprendre la machine $\hat{f}$, une part mise de côté, l’échantillon de test, de taille $m$, pour l’évaluer. Le risque de test est

$$\hat{R}_{\text{test}}(\hat{f}) = \frac{1}{m} \sum_{j=1}^{m} \ell\big(Y_j', \hat{f}(X_j')\big)$$

où les $(X_j', Y_j')$ sont les exemples de test, que la machine n’a jamais vus.

Pourquoi cette estimation est-elle honnête, alors que le risque d’apprentissage ne l’était pas ? Rappelez-vous l’exercice 3 du chapitre 2 : pour une fonction fixée à l’avance, le risque empirique est un estimateur sans biais du vrai risque. Or, vue de l’échantillon de test, $\hat{f}$ est justement fixée : elle a été choisie sans regarder ces exemples. Donc $\hat{R}_{\text{test}}(\hat{f})$ estime sans biais $R(\hat{f})$.

Sans biais ne veut pas dire exact. L’estimation varie selon les exemples tombés dans le test, et d’autant plus que $m$ est petit. C’est ici qu’intervient la construction de la confiance (étape 8 de la roue).

Un intervalle de confiance pour le taux d’erreur

En classification avec la perte 0-1, chaque exemple de test est une erreur ou non, avec une probabilité $R(\hat{f})$ : le nombre d’erreurs suit une loi binomiale. Pour $m$ assez grand, un intervalle de confiance à 95 % pour le vrai taux d’erreur est

$$\hat{e} \pm 1{,}96 \sqrt{\frac{\hat{e}\,(1 - \hat{e})}{m}}$$

• $\hat{e}$ est le taux d’erreur mesuré sur le test ;
• $m$ est le nombre d’exemples de test ;
• $1{,}96$ est le quantile de la loi normale qui laisse 2,5 % de chaque côté.

Exemple

Reprenons la classification « salaire supérieur à 60 k$ ». Sur 200 diplômés de test, le modèle se trompe 30 fois : $\hat{e} = 0{,}15$. La demi-largeur vaut $1{,}96 \sqrt{0{,}15 \times 0{,}85 / 200} \approx 0{,}049$. Le vrai taux d’erreur est donc, avec une confiance de 95 %, entre 0,10 et 0,20 environ. Pour gagner en précision, il faut beaucoup plus d’exemples : la largeur ne diminue qu’en $1/\sqrt{m}$, donc il faut quatre fois plus de test pour diviser l’incertitude par deux.

En régression, on procède de même avec la moyenne des pertes : si $s$ est l’écart-type des pertes individuelles $\ell(Y_j', \hat{f}(X_j'))$ sur le test, l’intervalle approché est $\hat{R}_{\text{test}} \pm 1{,}96\, s / \sqrt{m}$. Les bases de ces intervalles sont présentées dans le cours de Statistiques.

La largeur de l'intervalle décroît en 1/√m : passer de 200 à 800 exemples de test ne fait que diviser l'incertitude par deux.

À retenir

Un score sans intervalle de confiance est un demi-résultat. Avec 200 exemples de test, un taux d’erreur de 15 % signifie « entre 10 et 20 % » : deux modèles à 15 % et 17 % ne sont pas distinguables sur un tel échantillon.

La validation croisée

Mettre des données de côté a un coût : on apprend avec moins d’exemples, et l’estimation dépend du découpage. La validation croisée à $K$ plis (K-fold cross-validation) utilise chaque exemple tour à tour pour apprendre et pour évaluer : on coupe les données en $K$ blocs, on apprend sur $K - 1$ blocs, on évalue sur le bloc restant, et on moyenne les $K$ risques obtenus. Les détails pratiques (stratification, découpage par groupe ou dans le temps, pipelines contre les fuites de données) sont traités dans le cours d’Apprentissage automatique. Voyons plutôt ce que la validation croisée estime réellement.

Chaque bloc sert une fois de test et quatre fois à l'apprentissage. Le score de validation croisée est la moyenne des cinq risques mesurés sur des données non vues.

Chaque machine de la validation croisée est apprise sur une fraction $\frac{K-1}{K}$ des données, puis évaluée sur des exemples qu’elle n’a pas vus. Le score moyen estime donc le risque de la procédure « apprendre dans $\mathcal{H}$ avec environ $\frac{K-1}{K}\,n$ exemples ». Deux conséquences :

• l’estimation est légèrement pessimiste pour la machine finale, apprise sur les $n$ exemples (plus de données, moins de variance) ;
• elle évalue une méthode, pas une machine précise : les $K$ machines intermédiaires sont jetées après usage.

C’est l’outil idéal pour régler la capacité du chapitre 4, puisque le risque empirique d’apprentissage en est incapable. Sur nos diplômés, on estime par validation croisée le risque de chaque degré de polynôme et l’on garde le meilleur :

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(0)
n = 60
gpa = rng.uniform(2, 4, n)
salaire = 70 - 25 * np.exp(-2 * (gpa - 2)) + rng.normal(0, 4, n)   # population simulée

K = 5
plis = np.array_split(rng.permutation(n), K)          # K blocs d'indices, tirés au hasard

for degre in range(0, 7):
    risques = []
    for k in range(K):
        test = plis[k]
        appr = np.concatenate([plis[j] for j in range(K) if j != k])
        modele = Polynomial.fit(gpa[appr], salaire[appr], degre)      # apprendre sans le pli k
        risques.append(np.mean((salaire[test] - modele(gpa[test])) ** 2))  # évaluer sur le pli k
    print(f"degré {degre} : risque estimé {np.mean(risques):6.2f} (écart-type {np.std(risques):5.2f})")

Vous devriez voir le risque estimé baisser nettement du degré 0 au degré 2, puis stagner ou remonter : c’est le U du chapitre 4, cette fois reconstruit à partir d’un seul échantillon, sans connaître $f^*$.

Le résultat du code ci-dessus : le risque chute du degré 0 au degré 2, puis forme un palier proche du bruit σ² = 16. Au sein de ce palier, les écarts sont bien plus petits que la dispersion entre plis.

Piège

Choisir le meilleur degré par validation croisée, puis annoncer ce même score comme performance attendue. Ce score a servi à choisir : il est de nouveau optimiste (on a retenu le degré qui avait, entre autres, de la chance sur ces plis). Pour une estimation honnête après sélection, il faut soit un échantillon de test mis de côté dès le départ, soit une validation croisée imbriquée (nested) : une boucle interne choisit le degré, une boucle externe évalue la procédure complète.

Le bootstrap

Le bootstrap part d’une idée audacieuse : faire comme si l’échantillon était la population. On tire, avec remise, $n$ exemples parmi les $n$ disponibles. Certains diplômés apparaissent deux ou trois fois, d’autres pas du tout. On répète l’opération des centaines de fois, et chaque tirage joue le rôle d’un « nouvel échantillon ».

Combien d’exemples restent en dehors d’un tirage ? La probabilité qu’un diplômé donné ne soit jamais tiré en $n$ tirages est

$$\Big(1 - \frac{1}{n}\Big)^n \xrightarrow[n \to \infty]{} e^{-1} \approx 0{,}368$$

$n$	10	100	1 000	limite
$(1 - 1/n)^n$	0,349	0,366	0,368	0,368

Environ 36,8 % des exemples restent donc hors du tirage : on les appelle les exemples hors sac (out-of-bag). Ils fournissent un échantillon de test gratuit pour la machine apprise sur ce tirage. En moyennant les risques hors sac sur de nombreux tirages, on obtient une estimation du vrai risque, proche en esprit de la validation croisée. C’est d’ailleurs ce mécanisme qu’utilisent les forêts aléatoires pour s’auto-évaluer.

Trois tirages avec remise parmi 10 diplômés : les doublons prennent la place d'autres diplômés, qui restent hors sac et servent de test. Quand n grandit, la part hors sac tend vers e⁻¹ ≈ 36,8 %.

Le bootstrap sert aussi, et surtout, à mesurer l’incertitude de n’importe quelle statistique, même quand aucune formule n’existe. Par exemple, pour la pente $\hat{\beta}_1$ de la droite « salaire en fonction du GPA » :

import numpy as np

def pente(x, y):
    xc = x - x.mean()
    return (xc * (y - y.mean())).sum() / (xc ** 2).sum()

def intervalle_bootstrap(x, y, B=2000, seed=0):
    rng = np.random.default_rng(seed)
    n = len(x)
    pentes = []
    for _ in range(B):
        idx = rng.integers(0, n, n)          # n tirages avec remise
        pentes.append(pente(x[idx], y[idx]))
    return np.percentile(pentes, [2.5, 97.5])  # intervalle à 95 % par les percentiles

Si l’intervalle obtenu va de 5 à 17 k$ par point de GPA, vous savez que le « 11 k$ » annoncé au chapitre 2 est très incertain, ce que le chiffre seul ne disait pas.

Attention

Le bootstrap suppose que les exemples sont indépendants et tirés de la même loi. Si plusieurs lignes concernent un même diplômé (un salaire par année, par exemple), il faut tirer les diplômés avec remise, et non les lignes, sinon l’incertitude est sous-estimée.

Comparer au moins deux approches

L’étape 7 de la roue posait une règle : pour les mêmes données, au moins deux approches différentes, jugées sur leur performance sur la population plutôt que sur l’échantillon. Deux analystes n’interprètent pas le problème de la même façon, et c’est une richesse : la droite du premier et la méthode des plus proches voisins du second se contrôlent mutuellement.

Pour que la comparaison soit juste, il faut évaluer les deux approches sur les mêmes plis. On peut alors raisonner sur les différences pli par pli, ce qui élimine une grande partie du hasard lié au découpage (un pli « difficile » l’est pour les deux méthodes).

Exemple

Le service carrières compare, en classification « plus de 60 k$ », une régression logistique sur le GPA (A) et une méthode des plus proches voisins (B), sur les mêmes 5 plis (chiffres d’illustration) :

Pli	1	2	3	4	5	Moyenne
Erreur de A	0,18	0,21	0,17	0,20	0,19	0,190
Erreur de B	0,16	0,20	0,15	0,19	0,16	0,172
Différence A − B	0,02	0,01	0,02	0,01	0,03	0,018

Les erreurs de chaque méthode varient de 0,15 à 0,21 selon le pli, mais B fait mieux que A sur les cinq plis. La différence moyenne vaut 0,018, avec un écart-type de 0,0084 ; le rapport $0{,}018 / (0{,}0084 / \sqrt{5}) \approx 4{,}8$ est grand, ce qui plaide pour un avantage réel de B.

D'un pli à l'autre, les erreurs varient beaucoup (de 0,15 à 0,21), mais chaque segment descend : sur les mêmes plis, B bat A à chaque fois. Raisonner sur les différences élimine la difficulté propre à chaque pli.

Attention

Ce rapport ressemble à une statistique de Student appariée, mais les plis ne sont pas indépendants (les ensembles d’apprentissage se recouvrent largement). Le test habituel est donc trop optimiste. Pour une décision importante, on répète la validation croisée avec plusieurs découpages, on utilise un test corrigé pour ce recouvrement, ou l’on compare les deux modèles sur un vrai échantillon de test commun (avec un test de McNemar sur les exemples où ils sont en désaccord).

Et si aucune des deux approches ne se détache ? C’est aussi un résultat : on choisit alors la plus simple, la plus lisible ou la moins coûteuse à déployer. Au service carrières, une droite que chacun comprend vaut mieux qu’une méthode opaque qui ne gagne qu’un point d’erreur dans la marge d’incertitude.

En entretien

« Je ne me contente jamais d’un seul modèle : je compare au moins deux approches, dont une base simple, sur les mêmes plis de validation croisée, en regardant les différences pli par pli et leur variabilité. Et je garde un échantillon de test, utilisé une seule fois, pour annoncer la performance finale avec un intervalle de confiance. »

En résumé

• Le vrai risque étant incalculable, on l’estime par rééchantillonnage : échantillon de test, validation croisée, bootstrap.
• Le risque sur un échantillon de test jamais vu estime sans biais $R(\hat{f})$ ; on l’accompagne d’un intervalle de confiance, $\hat{e} \pm 1{,}96\sqrt{\hat{e}(1 - \hat{e})/m}$ pour un taux d’erreur.
• La validation croisée estime le risque d’une procédure d’apprentissage ; elle sert à régler la capacité, mais le score du meilleur réglage est optimiste (validation croisée imbriquée ou test séparé).
• Le bootstrap tire $n$ exemples avec remise ; environ 36,8 % restent hors du tirage et servent de test. Il mesure aussi l’incertitude de n’importe quelle statistique.
• On compare au moins deux approches sur les mêmes plis, en raisonnant sur les différences, sans oublier que les plis ne sont pas indépendants.

Exercices

Exercice · Un intervalle de confiance

facile

Un modèle se trompe sur 12 des 80 diplômés d’un échantillon de test. Donnez un intervalle de confiance approché à 95 % pour son vrai taux d’erreur. Combien faudrait-il d’exemples de test, à taux d’erreur égal, pour obtenir une demi-largeur d’environ 0,025 ?

Voir la solution

$\hat{e} = 12/80 = 0{,}15$. La demi-largeur vaut $1{,}96\sqrt{0{,}15 \times 0{,}85 / 80} \approx 0{,}078$. L’intervalle est donc environ $[0{,}07 \,;\, 0{,}23]$ : très large.

La demi-largeur est proportionnelle à $1/\sqrt{m}$. Pour passer de 0,078 à 0,025, il faut la diviser par un peu plus de 3, donc multiplier $m$ par environ 10 : avec $m = 800$, on trouve $1{,}96\sqrt{0{,}15 \times 0{,}85 / 800} \approx 0{,}025$. Il faut donc environ 800 exemples de test.

En code

Python

import numpy as np

erreurs, m = 12, 80
e = erreurs / m
demi = 1.96 * np.sqrt(e * (1 - e) / m)
print(e, round(demi, 3), np.round([e - demi, e + demi], 2))   # 0.15  0.078  [0.07 0.23]

# Taille de test pour une demi-largeur de 0,025 : m = 1,96² e(1 - e) / 0,025²
m_requis = 1.96 ** 2 * e * (1 - e) / 0.025 ** 2
print(round(m_requis))                            # 784, soit environ 800
print(round(1.96 * np.sqrt(e * (1 - e) / 800), 4))  # 0.0247

R

erreurs <- 12; m <- 80
e <- erreurs / m
demi <- 1.96 * sqrt(e * (1 - e) / m)
c(e, demi, e - demi, e + demi)                    # 0.15  0.078  0.072  0.228

# Taille de test pour une demi-largeur de 0,025
m_requis <- 1.96^2 * e * (1 - e) / 0.025^2
m_requis                                          # 783.6, soit environ 800
1.96 * sqrt(e * (1 - e) / 800)                    # 0.0247

Exercice · Compter les exemples hors sac

moyen

1. Avec $n = 10$ diplômés, quelle est la probabilité qu’un diplômé donné ne soit pas tiré dans un échantillon bootstrap ? Combien de diplômés, en moyenne, sont hors sac ?
2. Même question avec $n = 1\,000$.
3. Pourquoi l’estimation hors sac du risque est-elle, comme la validation croisée, un peu pessimiste ?

Voir la solution

1. $(1 - 1/10)^{10} = 0{,}9^{10} \approx 0{,}349$. En moyenne, $10 \times 0{,}349 \approx 3{,}5$ diplômés sont hors sac.
2. $(1 - 1/1000)^{1000} \approx 0{,}368$, soit environ 368 diplômés hors sac en moyenne.
3. Chaque machine bootstrap est apprise sur un tirage qui ne contient qu’environ 63,2 % d’exemples distincts (les autres sont des doublons). Elle dispose donc de moins d’information que la machine finale apprise sur tous les exemples, et son risque est en moyenne un peu plus élevé. L’estimation hors sac mesure le risque de ces machines « appauvries », donc elle surestime légèrement celui de la machine finale.

En code

Python

import numpy as np

for n in (10, 1000):
    p = (1 - 1 / n) ** n                          # P(un diplômé donné n'est jamais tiré)
    print(n, round(p, 3), round(n * p, 1))        # 10 : 0.349 3.5 / 1000 : 0.368 367.7

# Simulation : proportion d'exemples hors sac et d'exemples distincts
rng = np.random.default_rng(0)
n = 1000
hors_sac = [n - len(np.unique(rng.integers(0, n, n))) for _ in range(2000)]
print(np.mean(hors_sac) / n)                      # 0.368, donc 63,2 % de distincts

R

for (n in c(10, 1000)) {
  p <- (1 - 1 / n)^n                              # P(un diplômé donné n'est jamais tiré)
  print(c(n = n, p = p, hors_sac = n * p))        # 10 : 0.349 3.49 / 1000 : 0.368 367.7
}

# Simulation : proportion d'exemples hors sac et d'exemples distincts
set.seed(0)
n <- 1000
hors_sac <- replicate(2000, n - length(unique(sample.int(n, n, replace = TRUE))))
mean(hors_sac) / n                                # 0.367, donc 63,2 % de distincts

Exercice · Critiquer une comparaison

difficile

Un analyste compare une droite et un polynôme de degré réglable pour prédire les salaires. Il procède ainsi :

1. il évalue la droite par validation croisée à 5 plis, avec un premier découpage aléatoire ;
2. pour le polynôme, il essaie les degrés 1 à 10 par validation croisée à 5 plis, avec un second découpage aléatoire, et retient le meilleur degré ;
3. il compare le score de la droite au meilleur score du polynôme et conclut que le polynôme est meilleur.

Relevez les défauts de ce protocole et proposez-en un meilleur.

Voir la solution

• Des plis différents : les deux approches ne sont pas évaluées sur les mêmes découpages, donc une partie de l’écart peut venir du hasard du découpage. Il faut utiliser les mêmes plis pour les deux, et raisonner sur les différences pli par pli.
• Un score optimiste pour le polynôme : le score retenu est le meilleur de dix essais, choisi sur ces mêmes plis. La droite, elle, n’a eu qu’un essai. La comparaison favorise mécaniquement le polynôme. Il faut une validation croisée imbriquée (le choix du degré se fait à l’intérieur de chaque pli externe) ou un échantillon de test mis de côté au départ.
• Pas de mesure d’incertitude : sans écart-type ni intervalle, on ne sait pas si l’écart dépasse la variabilité naturelle.
• Détail révélateur : le degré 1 est la droite. Si le meilleur degré est 1, le polynôme ne peut pas être « meilleur » que la droite ; si c’est un autre degré, il faut vérifier que le gain est net et stable.

Protocole proposé : fixer un découpage externe en 5 plis ; dans chaque pli externe, choisir le degré par validation croisée interne sur la part d’apprentissage, puis évaluer la droite et le polynôme retenu sur le pli externe ; comparer les deux séries de cinq scores par leurs différences, et répéter l’ensemble avec plusieurs découpages.

En code

Python

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(0)
f_star = lambda x: 70 - 25 * np.exp(-2 * (x - 2))
def tirer(n):
    x = rng.uniform(2, 4, n)
    return x, f_star(x) + rng.normal(0, 4, n)
x_pop, y_pop = tirer(100_000)                     # sert à calculer le vrai risque

def score_cv(x, y, degre, K=5):
    plis = np.array_split(np.arange(len(x)), K)   # les x sont déjà dans un ordre aléatoire
    risques = []
    for k in range(K):
        appr = np.concatenate(plis[:k] + plis[k + 1:])
        m = Polynomial.fit(x[appr], y[appr], degre)
        risques.append(np.mean((y[plis[k]] - m(x[plis[k]])) ** 2))
    return np.mean(risques)

annonce, vrai = [], []
for _ in range(300):                              # 300 analystes, 300 échantillons de 60
    x, y = tirer(60)
    scores = [score_cv(x, y, d) for d in range(1, 11)]
    d = 1 + int(np.argmin(scores))                # le « meilleur » degré
    annonce.append(min(scores))                   # score annoncé par l'analyste
    m = Polynomial.fit(x, y, d)
    vrai.append(np.mean((y_pop - m(x_pop)) ** 2)) # vrai risque du polynôme retenu
print(round(np.mean(annonce), 2), round(np.mean(vrai), 2))   # 17.21  19.36 : optimiste

R

set.seed(0)
f_star <- function(x) 70 - 25 * exp(-2 * (x - 2))
tirer <- function(n) {
  x <- runif(n, 2, 4)
  data.frame(x = x, y = f_star(x) + rnorm(n, 0, 4))
}
pop <- tirer(100000)                              # sert à calculer le vrai risque

score_cv <- function(d, degre, K = 5) {
  pli <- rep(1:K, length.out = nrow(d))           # les lignes sont déjà dans un ordre aléatoire
  mean(sapply(1:K, function(k) {
    m <- lm(y ~ poly(x, degre), data = d[pli != k, ])
    mean((d$y[pli == k] - predict(m, d[pli == k, ]))^2)
  }))
}

res <- replicate(300, {                           # 300 analystes, 300 échantillons de 60
  d <- tirer(60)
  scores <- sapply(1:10, function(deg) score_cv(d, deg))
  m <- lm(y ~ poly(x, which.min(scores)), data = d)   # le « meilleur » degré
  c(annonce = min(scores), vrai = mean((pop$y - predict(m, pop))^2))
})
rowMeans(res)                                     # annonce 17.26, vrai 18.02 : optimiste