14 - Bootstrap e resampling

Hai quindici tempi di risposta di un servizio e vuoi dire quanto è incerta la loro mediana. Per la media esiste una formula, l'errore standard s/√n; per la mediana, per un rapporto fra due medie o per un indice di concentrazione la formula o non c'è o richiede ipotesi sulla distribuzione che non sai giustificare. Il bootstrap, proposto da Efron nel 1979, aggira il problema con un'idea concreta: se il campione è la miglior descrizione che hai della popolazione, estrai da esso nuovi campioni e guarda quanto varia la statistica. Questo articolo, il quattordicesimo della serie Probabilità e Statistica, costruisce il bootstrap a mano, confronta tre tipi di intervallo misurandone la copertura, mostra un caso in cui il metodo fallisce, e arriva al test di permutazione e al bootstrap a blocchi per dati dipendenti. Presuppone la distribuzione campionaria di uno stimatore, vista nell'articolo 6.

L'algoritmo: ricampionare con reinserimento

Dato un campione di n valori, un campione bootstrap si ottiene estraendo n valori dal campione con reinserimento: alcuni valori compaiono più volte, altri mai. Su ciascuno si ricalcola la statistica di interesse (qui la mediana), e ripetendo B volte si ottiene la distribuzione bootstrap, che approssima la distribuzione campionaria della statistica. La sua deviazione standard è l'errore standard bootstrap. Il listato genera quindici tempi di risposta da una distribuzione lognormale (asimmetrica, come sono spesso le latenze) con un seme fisso, e usa B = 9999. Da quella distribuzione ricava tre intervalli di confidenza al 95%. Il percentile prende i quantili 2,5% e 97,5% così come sono. Il basic li riflette attorno alla stima: 2t − q_97,5 e 2t − q_2,5. Il BCa (bias-corrected and accelerated, Efron 1987) corregge i quantili per la distorsione e per l'asimmetria dell'errore, e qui è calcolato da scipy.stats.bootstrap.

import numpy as np
from scipy import stats

rng = np.random.default_rng(2026)
# 15 tempi di risposta (ms) di un servizio: distribuzione asimmetrica, con code lunghe
x = np.round(rng.lognormal(mean=5.0, sigma=0.6, size=15)).astype(int)
print("dati ordinati:", np.sort(x))
n, B = len(x), 9999
med = np.median(x)

# Ricampiono n valori CON reinserimento, B volte, e ricalcolo la mediana su ciascun campione
idx = rng.integers(0, n, size=(B, n))
boot = np.median(x[idx], axis=1)
print(f"mediana osservata = {med:.1f}   errore standard bootstrap = {boot.std(ddof=1):.2f}")
print("quantili della distribuzione bootstrap (2.5%, 50%, 97.5%):", np.round(np.quantile(boot, [.025, .5, .975]), 1))

q_lo, q_hi = np.quantile(boot, [.025, .975])
perc = (q_lo, q_hi)                               # percentile: i quantili cosi' come sono
basic = (2 * med - q_hi, 2 * med - q_lo)          # basic: riflette i quantili attorno alla stima
bca = stats.bootstrap((x,), np.median, confidence_level=0.95, n_resamples=B, method="BCa", random_state=1).confidence_interval
print(f"percentile 95% [{perc[0]:.1f}, {perc[1]:.1f}]   basic [{basic[0]:.1f}, {basic[1]:.1f}]   BCa [{bca.low:.1f}, {bca.high:.1f}]")
print(f"quante mediane bootstrap sono valori distinti: {len(np.unique(boot))} su {B}")

# Quanto rumore aggiunge un B piccolo? Ripeto la stima dell'errore standard 200 volte per ogni B
for b in (50, 200, 1000, 5000):
    stime = [np.median(x[rng.integers(0, n, size=(b, n))], axis=1).std(ddof=1) for _ in range(200)]
    print(f"B = {b:5d}: errore standard stimato {np.mean(stime):6.2f}  (variabilita' fra le ripetizioni: sd {np.std(stime):.2f})")

# Output:
# dati ordinati: [ 48  92 123 125 126 130 141 143 163 171 178 202 218 229 343]
# mediana osservata = 143.0   errore standard bootstrap = 18.51
# quantili della distribuzione bootstrap (2.5%, 50%, 97.5%): [125. 143. 178.]
# percentile 95% [125.0, 178.0]   basic [108.0, 161.0]   BCa [126.0, 202.0]
# quante mediane bootstrap sono valori distinti: 13 su 9999
# B =    50: errore standard stimato  18.40  (variabilita' fra le ripetizioni: sd 1.94)
# B =   200: errore standard stimato  18.42  (variabilita' fra le ripetizioni: sd 0.99)
# B =  1000: errore standard stimato  18.46  (variabilita' fra le ripetizioni: sd 0.41)
# B =  5000: errore standard stimato  18.46  (variabilita' fra le ripetizioni: sd 0.20)

La mediana osservata è 143 ms e l'errore standard bootstrap 18,51. I tre intervalli sono diversi: percentile [125, 178], basic [108, 161], BCa [126, 202]. Non è un errore: sono costruzioni diverse che coincidono solo quando la distribuzione bootstrap è simmetrica e centrata sulla stima, e qui non lo è. C'è anche un segnale di allarme: le 9999 mediane bootstrap assumono solo 13 valori distinti, perché la mediana di 15 numeri può essere soltanto uno dei 15 valori osservati; con campioni così piccoli, la distribuzione bootstrap della mediana è «a gradini», e i quantili cadono su quei gradini. Nell'ultima parte il listato ripete la stima dell'errore standard 200 volte per ogni B: con B = 50 la variabilità fra ripetizioni è 1,94, con B = 1000 scende a 0,41 e con B = 5000 a 0,20, mentre la stima media (18,4-18,5) quasi non cambia. Il numero di ricampionamenti non rende la stima migliore, solo meno rumorosa: oltre qualche migliaio si paga tempo per ridurre un errore che è già piccolo rispetto a quello dovuto alle sole quindici osservazioni.

Quale intervallo? Misurare la copertura

Un intervallo al 95% dovrebbe contenere il valore vero nel 95% dei campioni. Con dati reali non lo sai mai, ma in simulazione sì: il programma qui sotto usa una popolazione lognormale con mediana vera 1 e media vera e^0,5 = 1,6487, estrae 1000 campioni di 30 osservazioni, costruisce gli intervalli per la mediana e per la media con ciascun metodo (B = 1000) e conta quante volte contengono il valore vero. Include anche l'intervallo t classico per la media come riferimento.

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
n, B, R = 30, 1000, 1000                       # taglia del campione, ricampionamenti, ripetizioni
vero = {"mediana": 1.0, "media": np.exp(0.5)}  # popolazione lognormale(0, 1): valori veri noti
statistiche = {"mediana": np.median, "media": np.mean}

copertura = {}
def conta(chiave, lo, hi, verita):
    copertura[chiave] = copertura.get(chiave, 0) + (lo <= verita <= hi)

for _ in range(R):
    x = rng.lognormal(0.0, 1.0, n)
    for nome, f in statistiche.items():
        t0 = f(x)
        boot = f(x[rng.integers(0, n, size=(B, n))], axis=1)   # B ricampionamenti con reinserimento
        lo, hi = np.quantile(boot, [.025, .975])
        conta((nome, "percentile"), lo, hi, vero[nome])
        conta((nome, "basic"), 2 * t0 - hi, 2 * t0 - lo, vero[nome])
        ci = stats.bootstrap((x,), f, n_resamples=B, method="BCa", random_state=rng).confidence_interval
        conta((nome, "BCa"), ci.low, ci.high, vero[nome])
    lo, hi = stats.t.interval(0.95, n - 1, loc=x.mean(), scale=stats.sem(x))   # intervallo t classico
    conta(("media", "t"), lo, hi, vero["media"])

print(f"copertura su {R} campioni (n={n}, B={B}); nominale 0.95, errore ~0.007")
for (s, m), c in copertura.items():
    print(f"  {s:8s} {m:10s} {c / R:.3f}")

# Output:
# copertura su 1000 campioni (n=30, B=1000); nominale 0.95, errore ~0.007
#   mediana  percentile 0.952
#   mediana  basic      0.814
#   mediana  BCa        0.954
#   media    percentile 0.871
#   media    basic      0.845
#   media    BCa        0.890
#   media    t          0.876

Per la mediana, percentile (0,952) e BCa (0,954) sono vicini al nominale 0,95 (l'errore di simulazione è circa 0,007), mentre il basic è a 0,814: il metodo che riflette i quantili può sbagliare di molto quando la distribuzione è asimmetrica. Per la media di dati lognormali con n = 30 nessun metodo arriva al 95%: BCa 0,890, t 0,876, percentile 0,871, basic 0,845. È il limite di tutti i metodi sulla media di dati molto asimmetrici con pochi dati: la distribuzione campionaria della media non è ancora vicina alla normale (articolo 5), e il bootstrap, che stima la distribuzione dai dati, eredita la stessa lentezza. In questo esperimento BCa è il migliore o il pari migliore in tutti i casi, ma non salva un campione troppo piccolo. Il bootstrap non crea informazione che i dati non contengono.

Dove il bootstrap fallisce: il massimo

Il bootstrap ordinario non funziona per statistiche che dipendono dagli estremi del campione. Il caso classico è il massimo di un campione uniforme su (0, θ), usato per stimare θ. In un ricampionamento da n valori, il massimo originale è presente con probabilità 1 − (1 − 1/n)ⁿ, che per n = 50 vale 0,6358 e tende a 1 − 1/e ≈ 0,632 per n grande: nel 63% dei ricampionamenti il massimo bootstrap è esattamente quello osservato. Nel mondo vero, invece, il massimo di un secondo campione indipendente coincide con quello del primo con probabilità zero. La distribuzione bootstrap ha un grosso atomo dove la distribuzione vera è continua, e non la approssima.

import numpy as np

rng = np.random.default_rng(11)

# 1) Dove il bootstrap non funziona: il massimo di un campione uniforme su (0, theta)
n = 50
x = rng.uniform(0, 10, n)
massimo = x.max()
boot_max = x[rng.integers(0, n, size=(100_000, n))].max(axis=1)
esatta = 1 - (1 - 1 / n) ** n                  # P(il massimo del campione e' nel ricampionamento)
print(f"P(massimo bootstrap = massimo osservato): esatta {esatta:.4f}   simulata {np.mean(boot_max == massimo):.4f}")
# nel mondo vero, il massimo di un altro campione e' quasi mai uguale a quello osservato
veri = rng.uniform(0, 10, size=(100_000, n)).max(axis=1)
print(f"nel mondo vero P(massimo = valore osservato) = {np.mean(veri == massimo):.4f}")

# Output:
# P(massimo bootstrap = massimo osservato): esatta 0.6358   simulata 0.6345
# nel mondo vero P(massimo = valore osservato) = 0.0000

Le due probabilità calcolate dal listato, quella esatta e la simulata, sono 0,6358 e 0,6345, e nel mondo vero la probabilità di coincidenza è 0,0000 su 100.000 campioni. La conseguenza pratica è che gli intervalli bootstrap per un massimo, un minimo o un quantile molto estremo non hanno la copertura dichiarata, qualunque sia B. Per queste statistiche servono altri approcci, come la distribuzione esatta del massimo o il bootstrap con sottocampioni di taglia m < n.

Il jackknife, il fratello più vecchio

Il jackknife (Quenouille, poi Tukey) ricalcola la statistica togliendo un dato alla volta, ottenendo n valori parziali. Da quei valori si stimano la distorsione e l'errore standard con formule chiuse, senza numeri casuali. Il listato verifica due proprietà esatte. La prima: applicato alla varianza con divisore n (distorta), la correzione jackknife restituisce esattamente la varianza con divisore n − 1: 1,749077 contro 1,749077. La seconda: l'errore standard jackknife della media coincide con s/√n, 0,341475 in entrambi i casi. Per statistiche lisce come queste il jackknife è esatto e veloce.

import numpy as np

rng = np.random.default_rng(12)

# 2) Jackknife: tolgo un dato alla volta
y = rng.lognormal(0, 1, 15)
m = len(y)
def jackknife(stat, y):
    parziali = np.array([stat(np.delete(y, i)) for i in range(len(y))])
    t0 = stat(y)
    bias = (len(y) - 1) * (parziali.mean() - t0)
    se = np.sqrt((len(y) - 1) / len(y) * np.sum((parziali - parziali.mean()) ** 2))
    return t0, bias, se

t0, bias, se = jackknife(lambda a: a.var(ddof=0), y)
print(f"varianza con 1/n: {t0:.6f}  bias jackknife {bias:.6f}  corretta {t0 - bias:.6f}  varianza con 1/(n-1): {y.var(ddof=1):.6f}")
print(f"errore standard della media: jackknife {jackknife(np.mean, y)[2]:.6f}   formula s/sqrt(n) {y.std(ddof=1) / np.sqrt(m):.6f}")

# 3) Errore standard della mediana (n = 15): la verita' e due stimatori, mediati su 300 campioni
vera_se = np.median(rng.lognormal(0, 1, size=(200_000, m)), axis=1).std()
jk, bs = [], []
for _ in range(300):
    z = rng.lognormal(0, 1, m)
    jk.append(jackknife(np.median, z)[2])
    bs.append(np.median(z[rng.integers(0, m, size=(500, m))], axis=1).std(ddof=1))
print(f"errore standard della mediana, vero {vera_se:.3f}")
print(f"  jackknife: media {np.mean(jk):.3f}, sd fra campioni {np.std(jk):.3f}, massimo {np.max(jk):.3f}")
print(f"  bootstrap: media {np.mean(bs):.3f}, sd fra campioni {np.std(bs):.3f}, massimo {np.max(bs):.3f}")

# Output:
# varianza con 1/n: 1.632471  bias jackknife -0.116605  corretta 1.749077  varianza con 1/(n-1): 1.749077
# errore standard della media: jackknife 0.341475   formula s/sqrt(n) 0.341475
# errore standard della mediana, vero 0.344
#   jackknife: media 0.326, sd fra campioni 0.265, massimo 1.947
#   bootstrap: media 0.399, sd fra campioni 0.191, massimo 1.137

Per la mediana, che non è una funzione liscia dei dati, il quadro cambia. L'errore standard vero (stimato con 200.000 campioni) è 0,344. Su 300 campioni di 15 osservazioni, il jackknife ha media 0,326 ma deviazione standard fra campioni 0,265 e un massimo di 1,947: è molto instabile. Il bootstrap ha media 0,399 (leggermente alta) ma più stabile, con deviazione 0,191. Sul jackknife per la mediana, quindi, non conviene fare affidamento, mentre il bootstrap è più regolare; nessuno dei due è preciso con 15 dati. Il titolo stesso dell'articolo di Efron (1979) presenta il bootstrap come una nuova lettura del jackknife.

Test di permutazione: il test quando l'ipotesi nulla è lo scambio

Il bootstrap stima una distribuzione campionaria; il test di permutazione costruisce direttamente la distribuzione di una statistica sotto l'ipotesi nulla. Se due metodi di imballaggio sono equivalenti, le etichette «A» e «B» attaccate ai 12 tempi sono arbitrarie, e ogni modo di distribuirle fra i dati è ugualmente probabile. Si ricalcola la differenza fra le medie per ciascuna redistribuzione e si vede quanto la differenza osservata è estrema. Con 6 e 6 pacchi le redistribuzioni sono 924 e si enumerano tutte, con 20 pacchi per gruppo sarebbero circa 138 miliardi e si campionano a caso, aggiungendo 1 a numeratore e denominatore perché un p-value di Monte Carlo non sia mai zero.

import itertools
import numpy as np

# Test di permutazione: due metodi di imballaggio, tempo (minuti) per 6 pacchi ciascuno
a = np.array([12.1, 14.3, 13.8, 15.2, 12.9, 14.7])
b = np.array([15.9, 14.8, 16.4, 13.9, 17.2, 15.5])
tutti = np.concatenate([a, b])
oss = a.mean() - b.mean()

# Esatto: tutte le C(12,6) = 924 maniere di scegliere quali 6 valori chiamare "A"
diff = np.array([tutti[list(i)].mean() - np.delete(tutti, list(i)).mean()
                 for i in itertools.combinations(range(12), 6)])
p_esatto = np.mean(np.abs(diff) >= abs(oss) - 1e-12)
print(f"differenza osservata {oss:.3f}   p esatto = {int(p_esatto * len(diff))}/{len(diff)} = {p_esatto:.4f}")

# Monte Carlo: 20000 permutazioni a caso (il +1 evita un p-value esattamente zero)
rng = np.random.default_rng(3)
mc = np.array([(lambda p: p[:6].mean() - p[6:].mean())(rng.permutation(tutti)) for _ in range(20_000)])
p_mc = (1 + np.sum(np.abs(mc) >= abs(oss) - 1e-12)) / (1 + len(mc))
print(f"p Monte Carlo = {p_mc:.4f} +- {np.sqrt(p_mc * (1 - p_mc) / len(mc)):.4f}")

# Output:
# differenza osservata -1.783   p esatto = 26/924 = 0.0281
# p Monte Carlo = 0.0269 +- 0.0011

La differenza osservata è −1,783 minuti. Solo 26 delle 924 permutazioni producono una differenza almeno così grande in valore assoluto: p esatto 26/924 = 0,0281. Il Monte Carlo con 20.000 permutazioni dà 0,0269 ± 0,0011, coerente. Il test è «esatto» nel senso che non usa approssimazioni, ma risponde a una domanda precisa: se i gruppi fossero scambiabili, quanto sarebbe sorprendente questo risultato? Non permette di dedurre nulla sulla popolazione oltre a quel confronto. La differenza fra permutazione e bootstrap per i test si riduce a questo: la permutazione ricampiona senza reinserimento e impone l'ipotesi nulla; il bootstrap ricampiona con reinserimento e, per un test, richiede di aver prima «spostato» i dati in modo che l'ipotesi nulla sia vera.

Quando i dati non sono indipendenti: il bootstrap a blocchi

Il bootstrap ordinario presuppone osservazioni indipendenti: mescolare i dati di una serie temporale distrugge l'ordine e quindi la dipendenza. Il bootstrap a blocchi (Künsch, 1989) ricampiona blocchi consecutivi di lunghezza L, che conservano la dipendenza al loro interno. Il listato usa una serie autoregressiva AR(1) con coefficiente ρ = 0,7, in cui ogni valore dipende dal precedente: il suo errore standard vero della media è noto per simulazione (20.000 serie) e anche per formula asintotica, √(σ_x²/n · (1 + ρ)/(1 − ρ)). Si confrontano poi gli errori standard stimati dal bootstrap per cinque lunghezze di blocco, L = 1 essendo il bootstrap ordinario.

import numpy as np

# Serie temporale AR(1): x_t = rho * x_{t-1} + rumore. Gli osservati NON sono indipendenti.
rho, n = 0.7, 200
rng = np.random.default_rng(21)

def serie(r, k, burn=200):
    e = r.standard_normal((k, n + burn))
    x = np.zeros((k, n + burn))
    for t in range(1, n + burn):
        x[:, t] = rho * x[:, t - 1] + e[:, t]
    return x[:, burn:]

# Verita': errore standard della media su 20.000 serie indipendenti, e valore teorico asintotico
medie = serie(rng, 20_000).mean(axis=1)
sigma2_x = 1 / (1 - rho**2)
teorico = np.sqrt(sigma2_x / n * (1 + rho) / (1 - rho))
print(f"errore standard della media: simulato {medie.std():.4f}   teorico asintotico {teorico:.4f}   se i dati fossero indipendenti {np.sqrt(sigma2_x / n):.4f}")

def se_bootstrap(x, L, B=500):
    k = int(np.ceil(n / L))                      # numero di blocchi consecutivi di lunghezza L
    partenze = rng.integers(0, n - L + 1, size=(B, k))
    campioni = np.concatenate([x[partenze + j][:, :, None] for j in range(L)], axis=2).reshape(B, -1)[:, :n]
    return campioni.mean(axis=1).std(ddof=1)

campioni = serie(rng, 200)
for L in (1, 5, 10, 20, 40):                      # L = 1 e' il bootstrap ordinario
    stime = [se_bootstrap(x, L) for x in campioni]
    print(f"blocco L = {L:2d}: errore standard stimato (media su 200 serie) {np.mean(stime):.4f}  sd {np.std(stime):.4f}")

# Output:
# errore standard della media: simulato 0.2319   teorico asintotico 0.2357   se i dati fossero indipendenti 0.0990
# blocco L =  1: errore standard stimato (media su 200 serie) 0.0980  sd 0.0091
# blocco L =  5: errore standard stimato (media su 200 serie) 0.1689  sd 0.0227
# blocco L = 10: errore standard stimato (media su 200 serie) 0.1932  sd 0.0334
# blocco L = 20: errore standard stimato (media su 200 serie) 0.2018  sd 0.0445
# blocco L = 40: errore standard stimato (media su 200 serie) 0.1886  sd 0.0591

L'errore standard vero è 0,2319, e la formula asintotica 0,2357. Se i dati fossero indipendenti sarebbe 0,0990: ignorare la dipendenza sottostima l'incertezza di 2,3 volte, e il bootstrap ordinario (L = 1) infatti dà 0,0980. Con blocchi crescenti la stima sale a 0,1689 (L = 5), 0,1932 (L = 10) e 0,2018 (L = 20), poi ridiscende a 0,1886 con L = 40, quando i blocchi sono pochi e la variabilità fra serie aumenta (sd 0,0591). Nessuna lunghezza raggiunge il valore vero: il bootstrap a blocchi con blocchi mobili tende a sottostimare con n = 200. Il punto non è la precisione, ma la direzione: la scelta di L è un compromesso fra conservare la dipendenza (blocchi lunghi) e avere abbastanza blocchi diversi (blocchi corti).

Esercizio svolto e prova da fare da solo

Quanto vale la probabilità che il massimo di un campione di n = 10 compaia in un ricampionamento bootstrap? Ogni estrazione sceglie il massimo con probabilità 1/10, quindi lo manca con probabilità 9/10; in 10 estrazioni indipendenti lo manca sempre con probabilità (9/10)¹⁰ = 0,3487, e lo include con probabilità 1 − 0,3487 = 0,6513. Verificato con il calcolo in Python: 1 - 0.9**10 restituisce 0,6513215599.

Che cosa sai fare adesso

Sai costruire un bootstrap a mano, scegliere fra intervallo percentile, basic e BCa sapendo che il basic può avere copertura molto bassa su dati asimmetrici, e riconoscere i casi in cui non funziona, come il massimo di un campione. Sai applicare jackknife, test di permutazione e bootstrap a blocchi, con le limitazioni misurate qui sopra. Il bootstrap è l'alternativa frequentista agli intervalli di credibilità dell'articolo 11; per i test senza ipotesi di distribuzione vedi anche i test non parametrici. L'ultimo articolo della serie, l'analisi esplorativa dei dati, è il passo che precede qualunque di questi calcoli.