05 - Cross-validation e selezione del modello

Hai addestrato un modello e vuoi sapere quanto sbaglierà su dati nuovi. Il modo più immediato è tenere da parte il 20% degli esempi e misurarci sopra l'errore; il difetto è che quel numero dipende da quali esempi sono capitati nel 20%. La validazione incrociata (cross-validation) corregge il difetto ripetendo la misura su più partizioni, ma ha le sue trappole: ordine dei dati, esempi imparentati, tempo, e la tentazione di usare lo stesso punteggio sia per scegliere il modello sia per dichiararne la qualità. Questo articolo mostra con esperimenti riproducibili dove ciascuna di queste trappole morde e quale divisione la evita. Parte dalle idee dei capitoli 01 (baseline, insieme di test, leakage) e 04 (bias e varianza di un modello).

Gli esempi girano con scikit-learn 1.9.0 e numpy 2.4.6 su dataset inclusi (Breast Cancer, 569 esempi; Iris, 150) o su dati sintetici con seme fisso. Ogni uscita è stata ottenuta due volte con lo stesso risultato.

Perché una sola partizione non basta

La K-fold cross-validation divide i dati in K blocchi (fold) di dimensione simile; a turno ogni blocco fa da test e gli altri K − 1 da training; la stima finale è la media dei K punteggi. Con K uguale al numero di esempi si ottiene la variante leave-one-out. Il manuale di James, Witten, Hastie e Tibshirani (An Introduction to Statistical Learning) descrive il metodo nel paragrafo 5.1.3 e discute la scelta di K nel 5.1.4: con K uguale a 5 o 10 i training set contengono circa (K − 1)/K dei dati, quindi la stima è appena pessimista, mentre il leave-one-out ha bias minimo ma varianza più alta, perché mediare modelli addestrati su insiemi quasi identici dà risultati molto correlati. Gli autori concludono che 5 o 10 sono scelte che empiricamente non soffrono né di bias eccessivo né di varianza elevata.

Il listato implementa il K-fold a mano, lo confronta con cross_val_score, e misura quanto varia una stima cambiando soltanto la partizione. Il modello è la regressione logistica standardizzata del capitolo 03, nella Pipeline che ricalcola la standardizzazione dentro ogni fold.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
modello = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000))

# K-fold a mano: ogni blocco fa da test una volta, gli altri K-1 da training
def kfold_a_mano(K, seme):
    idx = np.random.RandomState(seme).permutation(len(y))      # stessa permutazione di KFold(shuffle=True)
    blocchi = np.array_split(idx, K)
    acc = []
    for k in range(K):
        train = np.concatenate([b for j, b in enumerate(blocchi) if j != k])
        acc.append(modello.fit(X[train], y[train]).score(X[blocchi[k]], y[blocchi[k]]))
    return np.mean(acc)

a_mano = kfold_a_mano(5, 0)
sk = cross_val_score(modello, X, y, cv=KFold(5, shuffle=True, random_state=0)).mean()
print("5-fold a mano:", round(a_mano, 4), " cross_val_score:", round(sk, 4), " uguali:", np.isclose(a_mano, sk))

# quanto varia la stima cambiando solo la partizione?
hold = []
for s in range(100):
    Xa, Xb, ya, yb = train_test_split(X, y, test_size=0.2, random_state=s)
    hold.append(modello.fit(Xa, ya).score(Xb, yb))
cv5 = [cross_val_score(modello, X, y, cv=KFold(5, shuffle=True, random_state=s)).mean() for s in range(100)]
print(f"hold-out 80/20, 100 partizioni: min={min(hold):.3f} max={max(hold):.3f} dev.std={np.std(hold):.4f}")
print(f"5-fold, 100 partizioni:         min={min(cv5):.3f} max={max(cv5):.3f} dev.std={np.std(cv5):.4f}")

# Output:
# 5-fold a mano: 0.9719  cross_val_score: 0.9719  uguali: True
# hold-out 80/20, 100 partizioni: min=0.947 max=1.000 dev.std=0.0138
# 5-fold, 100 partizioni:         min=0.967 max=0.984 dev.std=0.0037

La versione a mano e quella di scikit-learn coincidono (0,9719) perché usano la stessa permutazione, con blocchi di dimensioni uguali a meno di uno. Il risultato che conta è l'ultimo: ripetendo 100 volte una divisione 80/20 l'accuratezza stimata va da 0,947 a 1,000, con deviazione standard 0,0138; ripetendo 100 volte una 5-fold, da 0,967 a 0,984, con deviazione standard 0,0037. Con una divisione sola potresti dichiarare «100%» o «94,7%» per lo stesso identico modello, a seconda del seme. La 5-fold usa tutti i 569 esempi sia in training sia in test a turno e mediare riduce la deviazione standard di circa 3,7 volte in questo esperimento; il costo è fare K adattamenti invece di uno.

Quando il K-fold ingenuo inganna

Il K-fold dà una stima onesta solo se ogni blocco assomiglia a «dati nuovi» tanto quanto i dati futuri. Tre situazioni violano l'ipotesi, e ciascuna ha uno schema di divisione pensato per lei.

Dati ordinati e classi rare: stratificare

Il dataset Iris ha le 150 righe ordinate per specie: le prime 50 sono tutte della specie 0, le ultime 50 della specie 2. Un K-fold a 3 blocchi senza mescolare prende ogni blocco di 50 righe come test: il blocco di test contiene una specie che il modello, addestrato sulle altre due, non ha mai visto.

import numpy as np
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, StratifiedKFold, cross_val_score

X, y = load_iris(return_X_y=True)                  # le 150 righe sono ordinate per specie
m = LogisticRegression(max_iter=1000)
print("etichette nelle prime/ultime righe:", y[:3].tolist(), y[-3:].tolist())
print("KFold(3) senza mescolare:   ", cross_val_score(m, X, y, cv=KFold(3)).round(3).tolist())
print("KFold(3) mescolando:        ", cross_val_score(m, X, y, cv=KFold(3, shuffle=True, random_state=0)).round(3).tolist())
print("StratifiedKFold(3):         ", cross_val_score(m, X, y, cv=StratifiedKFold(3)).round(3).tolist())
print("cv=3 (default per classificatori):", cross_val_score(m, X, y, cv=3).round(3).tolist())

# Output:
# etichette nelle prime/ultime righe: [0, 0, 0] [2, 2, 2]
# KFold(3) senza mescolare:    [0.0, 0.0, 0.0]
# KFold(3) mescolando:         [0.98, 0.96, 0.96]
# StratifiedKFold(3):          [0.98, 0.96, 0.98]
# cv=3 (default per classificatori): [0.98, 0.96, 0.98]

Il risultato è 0,0 su tutti e tre i fold: il modello sbaglia ogni esempio perché non conosce la classe che deve prevedere. Mescolando (shuffle=True) o usando StratifiedKFold, che mantiene in ogni fold la stessa proporzione di classi del dataset, l'accuratezza torna a 0,96-0,98. La documentazione di scikit-learn precisa che quando cv è un intero, cross_val_score usa già StratifiedKFold per i classificatori (ultima riga), mentre una KFold passata esplicitamente sovrascrive questa scelta: l'errore della prima riga nasce dalla KFold senza mescolare. La stratificazione conta soprattutto con classi rare, dove un fold potrebbe non contenerne nessun esempio.

Esempi imparentati: raggruppare

Se il dataset contiene più righe per la stessa persona, macchina o cliente, un K-fold casuale mette righe della stessa persona sia in training sia in test, e il modello può riconoscere la persona invece di imparare il fenomeno. Il listato costruisce il caso estremo: 40 pazienti con 10 visite ciascuno, un tratto personale che si ripete in ogni visita, e un esito scelto a caso per paziente, quindi senza alcun legame con le misure. Il valore onesto è 0,5.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GroupKFold, KFold, cross_val_score

rng = np.random.default_rng(0)
n_pazienti, n_visite = 40, 10
impronta = rng.standard_normal((n_pazienti, 8))                    # tratto proprio di ogni paziente
X = np.repeat(impronta, n_visite, axis=0) + 0.3 * rng.standard_normal((n_pazienti * n_visite, 8))
esito = rng.integers(0, 2, n_pazienti)                             # esito per paziente, scelto a caso: nessun legame con le misure
y = np.repeat(esito, n_visite)
gruppi = np.repeat(np.arange(n_pazienti), n_visite)

rf = RandomForestClassifier(n_estimators=100, random_state=0)
casuale = cross_val_score(rf, X, y, cv=KFold(5, shuffle=True, random_state=0)).mean()
per_paziente = cross_val_score(rf, X, y, cv=GroupKFold(5), groups=gruppi).mean()
print("400 visite di 40 pazienti, esito casuale (limite onesto: 0.5)")
print("KFold casuale:", round(casuale, 3))
print("GroupKFold per paziente:", round(per_paziente, 3))

# Output:
# 400 visite di 40 pazienti, esito casuale (limite onesto: 0.5)
# KFold casuale: 0.95
# GroupKFold per paziente: 0.508

Con il K-fold casuale l'accuratezza stimata è 0,95 su un problema in cui nulla si può prevedere: il modello ha riconosciuto i pazienti dal tratto personale e ne ha ricordato l'esito. Con GroupKFold, che tiene tutte le visite di un paziente nello stesso fold e passa l'identificativo con groups=, la stima torna a 0,508. La documentazione di scikit-learn descrive lo stesso scenario, soggetti con più campioni e modelli abbastanza flessibili da imparare tratti personali. Il difetto è lo stesso del leakage del capitolo 01: informazione che in produzione non avresti, perché i pazienti nuovi non hanno visite nel training.

Dati nel tempo: non guardare il futuro

Quando le righe sono ordinate nel tempo e vuoi prevedere il futuro, un blocco di test nel passato rispetto al training è un esame truccato. TimeSeriesSplit fa sì che ogni test venga sempre dopo il suo training.

import numpy as np
from sklearn.model_selection import TimeSeriesSplit

for k, (tr, te) in enumerate(TimeSeriesSplit(n_splits=4).split(np.arange(12)), 1):
    print(f"split {k}: train = {tr.tolist()}  test = {te.tolist()}")

# Output:
# split 1: train = [0, 1, 2, 3]  test = [4, 5]
# split 2: train = [0, 1, 2, 3, 4, 5]  test = [6, 7]
# split 3: train = [0, 1, 2, 3, 4, 5, 6, 7]  test = [8, 9]
# split 4: train = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]  test = [10, 11]

I training set crescono e contengono sempre tutto ciò che precede il test: nello split 3 il modello si addestra sulle righe da 0 a 7 e viene valutato su 8 e 9. Come dice la documentazione, a differenza degli altri schemi i training set successivi sono sovrainsiemi dei precedenti. Il prezzo è che i primi split addestrano su pochissimi dati, e che non si può mescolare: se i dati hanno una stagionalità o un cambio di regime, la stima è comunque quella del passato.

Scegliere gli iperparametri

Un iperparametro è una scelta che non si impara dai dati ma si fissa prima: per la regressione logistica la forza della penalità, C (più piccolo, più penalità e modello più semplice). Si prova una griglia di valori, si calcola per ciascuno il punteggio in validazione incrociata e si prende il migliore: è ciò che fa GridSearchCV, che alla fine riadatta (refit) il modello migliore su tutto il training. Prima della ricerca si mette da parte un insieme di test, che si usa una volta sola.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=0, stratify=y)
pipe = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000))
Cs = np.logspace(-3, 2, 11)
gs = GridSearchCV(pipe, {"logisticregression__C": Cs}, cv=StratifiedKFold(5, shuffle=True, random_state=0)).fit(Xtr, ytr)

media, dev = gs.cv_results_["mean_test_score"], gs.cv_results_["std_test_score"]
i = int(np.argmax(media))
se = dev[i] / np.sqrt(5)                               # errore standard della stima del migliore
# regola dell'errore standard: il modello piu' semplice (C piu' piccolo) entro 1 SE dal migliore
j = min(k for k in range(len(Cs)) if media[k] >= media[i] - se)
print(f"C migliore: {Cs[i]:.3g}  accuratezza CV = {media[i]:.4f} (SE {se:.4f})")
print(f"C con la regola 1-SE: {Cs[j]:.3g}  accuratezza CV = {media[j]:.4f}")
print(f"test, una volta sola, con C migliore: {gs.score(Xte, yte):.4f}")

# Output:
# C migliore: 0.316  accuratezza CV = 0.9812 (SE 0.0054)
# C con la regola 1-SE: 0.1  accuratezza CV = 0.9765
# test, una volta sola, con C migliore: 0.9650

Il valore migliore è C = 0,316 con accuratezza di validazione 0,9812, e il suo errore standard è 0,0054 (deviazione standard dei cinque punteggi divisa per la radice di 5). La regola dell'errore standard sceglie invece il modello più semplice entro un errore standard dal migliore: qui C = 0,1 con 0,9765, quasi indistinguibile dal massimo ma con più penalità. Il manuale citato la descrive nel paragrafo 6.1.3, con la motivazione che fra modelli più o meno equivalenti conviene il più semplice, e conviene quando la differenza fra due modelli è sotto la soglia del rumore della stima. L'accuratezza sull'insieme di test, calcolata una sola volta, è 0,9650: più bassa di quella di validazione, come ci si aspetta da un valore che è stato scelto proprio perché massimo sulla validazione, e su 143 esempi sarebbe un errore ritenerla diversa in modo significativo.

La validazione annidata

Il punteggio migliore di una ricerca a griglia è il massimo di molte stime rumorose, quindi è sistematicamente ottimista: la configurazione vincente ha ricevuto parte della sua fortuna dai fold in cui è stata valutata. Per misurare l'ottimismo costruiamo un caso in cui la risposta onesta è nota: 60 esempi e 40 feature di puro rumore, etichette senza legame con le misure. Ripetiamo su 10 dataset diversi, confrontando il punteggio migliore della ricerca con la validazione annidata: un ciclo esterno divide i dati, e dentro ogni training esterno una ricerca a griglia indipendente sceglie gli iperparametri; il fold esterno, mai visto dalla ricerca, dà la stima onesta.

import numpy as np
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold, cross_val_score
from sklearn.pipeline import Pipeline

griglia = {"sel__k": [1, 2, 3, 5, 8, 12, 20, 40], "clf__C": [0.01, 0.1, 1, 10, 100]}   # 40 configurazioni
interno = StratifiedKFold(5, shuffle=True, random_state=0)
esterno = StratifiedKFold(5, shuffle=True, random_state=1)
migliore, annidata = [], []
for seme in range(10):
    rng = np.random.default_rng(seme)
    X = rng.standard_normal((60, 40))                  # 60 esempi, 40 feature di puro rumore
    y = np.repeat([0, 1], 30)                          # etichette senza alcun legame con X: limite onesto = 0.5
    pipe = Pipeline([("sel", SelectKBest(f_classif)), ("clf", LogisticRegression(max_iter=1000))])
    gs = GridSearchCV(pipe, griglia, cv=interno).fit(X, y)
    migliore.append(gs.best_score_)
    annidata.append(cross_val_score(GridSearchCV(pipe, griglia, cv=interno), X, y, cv=esterno).mean())
print("media su 10 dataset di rumore (valore vero: 0.5)")
print("miglior punteggio della ricerca (best_score_):", round(np.mean(migliore), 3))
print("validazione annidata:                         ", round(np.mean(annidata), 3))

# Output:
# media su 10 dataset di rumore (valore vero: 0.5)
# miglior punteggio della ricerca (best_score_): 0.588
# validazione annidata:                          0.49

Sulla griglia di 40 configurazioni (8 valori di k per SelectKBest e 5 valori di C) il miglior punteggio medio è 0,588, ma il valore vero è 0,5 perché i dati sono rumore: quasi nove punti percentuali sopra il valore vero, che sembrano un segnale. La validazione annidata dà 0,490, cioè il caso. Il calcolo costa: con 5 fold esterni e una griglia di 11 valori valutata su 5 fold interni sono 280 adattamenti, contro 56 della ricerca semplice (conteggi nell'esercizio sotto). La conclusione operativa è di usare la validazione semplice per scegliere gli iperparametri, ma non riportare il suo punteggio come stima delle prestazioni: o si usa l'annidata, o si usa un insieme di test mai toccato.

Esercizio svolto: quanti adattamenti, e su quanti esempi

Con i 569 esempi del dataset, una 5-fold e una griglia di 11 valori di C, calcola le dimensioni dei fold e il numero di adattamenti della ricerca semplice e di quella annidata.

import numpy as np
from sklearn.model_selection import KFold

n, K = 569, 5
dimensioni = [(len(tr), len(te)) for tr, te in KFold(K).split(np.zeros(n))]
print("(train, test) per ogni fold:", dimensioni)
configurazioni = 11                                   # valori di C nella griglia
print("adattamenti di una ricerca a griglia con CV a 5 blocchi:", configurazioni * K, "+ 1 di refit =", configurazioni * K + 1)
print("adattamenti di una validazione annidata 5 esterni x ricerca:", K * (configurazioni * K + 1))

# Output:
# (train, test) per ogni fold: [(455, 114), (455, 114), (455, 114), (455, 114), (456, 113)]
# adattamenti di una ricerca a griglia con CV a 5 blocchi: 55 + 1 di refit = 56
# adattamenti di una validazione annidata 5 esterni x ricerca: 280

569 non è divisibile per 5: quattro fold hanno 114 esempi di test e uno ne ha 113, e di conseguenza i training set sono di 455 e 456 esempi. La ricerca semplice costa 11 × 5 = 55 adattamenti più uno di refit; la annidata ripete l'intera ricerca in ognuno dei 5 fold esterni, quindi 5 × 56 = 280. Quando il modello è lento, il fattore cinque della annidata decide se sia praticabile, e la risposta realistica è spesso un insieme di test unico e di dimensione adeguata.

Quando la validazione incrociata non basta

• I dati futuri avranno una distribuzione diversa da quelli raccolti: la stima riguarda il passato, e serve un monitoraggio (capitolo 12).
• Gli esempi non sono indipendenti in modi che non conosci: serve conoscere la struttura (gruppi, tempo) per scegliere lo schema.
• La metrica non rispecchia il costo degli errori: una stima precisa di un numero sbagliato resta sbagliata (capitolo 11).

Le tre voci hanno in comune che la validazione incrociata garantisce solo una stima dell'errore su dati simili a quelli già visti; la somiglianza è un'ipotesi del progettista, non una proprietà del metodo.

Che cosa sai fare adesso

Sai implementare un K-fold e spiegare perché una sola partizione dà stime che saltano di qualche punto, scegliere fra stratificato, per gruppi e per tempo guardando la struttura dei dati, scegliere un iperparametro con la ricerca a griglia e la regola dell'errore standard, e riconoscere quando il punteggio di una ricerca è gonfiato e serve la validazione annidata. Il passo successivo è sfruttare gli alberi, che hanno molti iperparametri da validare: 06 - Decision Trees e Random Forest.

Nella stessa serie trovi anche: 02 - Regressione Lineare from scratch in Python e 07 - Support Vector Machines: kernel trick e SVM non lineari.