Metodi bootstrap per la deviazione standard

Bootstrap: la rivoluzione statistica dell’era informatica

Il ricampionamento bootstrap è una tecnica statistica potente che stima la distribuzione campionaria di qualsiasi statistica ricampionando ripetutamente dai dati osservati. Introdotto da Bradley Efron nel 1979, ha rivoluzionato l’inferenza statistica consentendo l’analisi di statistiche complesse senza dipendere da formule matematiche o assunzioni distributive.

L’intuizione alla base del bootstrap è elegantemente semplice: il campione è la migliore stima della popolazione. Ricampionando dal proprio campione (con reinserimento), si simula ciò che accadrebbe se si potesse campionare ripetutamente dalla popolazione. Questo approccio è particolarmente prezioso per la deviazione standard, dove le formule tradizionali per gli intervalli di confidenza assumono la normalità, un’assunzione che spesso non regge nella pratica.

Il bootstrap è diventato essenziale nella data science moderna perché funziona con qualsiasi statistica (mediana, correlazione, coefficienti di regressione, pesi delle reti neurali) e non fa assunzioni sulla distribuzione sottostante dei dati.

Perché il bootstrap per la deviazione standard?

Gli intervalli di confidenza tradizionali per la deviazione standard assumono che i dati provengano da una distribuzione normale. Quando questa assunzione viene meno (come spesso accade), questi intervalli possono essere gravemente imprecisi. Il bootstrap fornisce un’alternativa libera da assunzioni distributive.

Vantaggi chiave del bootstrap per la deviazione standard:

• Nessuna assunzione distributiva: Funziona ugualmente bene con dati normali, asimmetrici o a code pesanti
• Prestazioni con piccoli campioni: Spesso più accurato dei metodi parametrici con n < 30
• Gestisce statistiche complesse: Lo stesso approccio funziona per DS troncata, MAD o misure di variabilità personalizzate
• Insight visivo: La distribuzione bootstrap mostra ciò che accade, non solo i numeri finali

La procedura bootstrap

L’algoritmo bootstrap è notevolmente semplice. Dal campione originale di n osservazioni:

Quanti campioni bootstrap? B = 1.000 è spesso sufficiente per stime approssimative e test di ipotesi. Per intervalli di confidenza, B = 10.000 fornisce percentili stabili. Per intervalli BCa di qualità pubblicabile, si raccomandano B = 15.000 o più.

Metodi per gli intervalli di confidenza bootstrap

Esistono diversi metodi per costruire intervalli di confidenza dai campioni bootstrap, ciascuno con i propri compromessi:

1. Metodo dei percentili (il più semplice)

L’approccio più intuitivo: prendere direttamente i percentili della distribuzione bootstrap.

Per 10.000 campioni bootstrap, sono il 250° e il 9.750° valore ordinato. Semplice ma può essere distorto quando la distribuzione bootstrap è asimmetrica.

2. Bootstrap di base (pivotale)

Utilizza la relazione tra la statistica campionaria e le statistiche bootstrap:

Dove θ̂ è la DS campionaria originale. Questo “riflette” l’intervallo dei percentili attorno alla stima campionaria.

3. BCa (corretto per distorsione e accelerato)

Lo standard di riferimento per l’accuratezza. Il BCa corregge sia la distorsione nella distribuzione bootstrap che l’accelerazione (come l’errore standard cambia con il valore del parametro). Più complesso da calcolare ma fornisce intervalli accurati al secondo ordine.

Esempio svolto: dati non normali

Consideriamo 15 misurazioni di tempi di risposta (in ms): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312. Questi dati sono asimmetrici a destra (alcune risposte molto lente).

L’IC bootstrap è asimmetrico (più ampio sul lato alto), riflettendo la natura asimmetrica a destra dei dati. L’IC chi-quadrato non cattura questa asimmetria.

Implementazione in Python

Implementazione completa del bootstrap con più metodi per gli IC:

import numpy as np
from scipy import stats

def bootstrap_sd_ci(data, n_bootstrap=10000, ci=0.95, method='percentile'):
    """
    Bootstrap confidence interval for standard deviation.

    Parameters:
    -----------
    data : array-like - Original sample
    n_bootstrap : int - Number of bootstrap samples
    ci : float - Confidence level (e.g., 0.95)
    method : str - 'percentile', 'basic', or 'bca'

    Returns:
    --------
    tuple : (lower_bound, upper_bound, bootstrap_sds)
    """
    data = np.array(data)
    n = len(data)
    original_sd = np.std(data, ddof=1)

    # Generate bootstrap samples and calculate SDs
    bootstrap_sds = np.array([
        np.std(np.random.choice(data, size=n, replace=True), ddof=1)
        for _ in range(n_bootstrap)
    ])

    alpha = 1 - ci

    if method == 'percentile':
        lower = np.percentile(bootstrap_sds, 100 * alpha/2)
        upper = np.percentile(bootstrap_sds, 100 * (1 - alpha/2))

    elif method == 'basic':
        lower = 2*original_sd - np.percentile(bootstrap_sds, 100*(1-alpha/2))
        upper = 2*original_sd - np.percentile(bootstrap_sds, 100*alpha/2)

    elif method == 'bca':
        # Bias correction
        prop_less = np.mean(bootstrap_sds < original_sd)
        z0 = stats.norm.ppf(prop_less)

        # Acceleration (jackknife estimate)
        jackknife_sds = np.array([
            np.std(np.delete(data, i), ddof=1) for i in range(n)
        ])
        jack_mean = jackknife_sds.mean()
        a = np.sum((jack_mean - jackknife_sds)**3) / \
            (6 * np.sum((jack_mean - jackknife_sds)**2)**1.5)

        # Adjusted percentiles
        z_alpha = stats.norm.ppf([alpha/2, 1-alpha/2])
        adj_percentiles = stats.norm.cdf(
            z0 + (z0 + z_alpha) / (1 - a*(z0 + z_alpha))
        ) * 100
        lower = np.percentile(bootstrap_sds, adj_percentiles[0])
        upper = np.percentile(bootstrap_sds, adj_percentiles[1])

    return lower, upper, bootstrap_sds

# Example usage
response_times = [245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312]

for method in ['percentile', 'basic', 'bca']:
    lower, upper, _ = bootstrap_sd_ci(response_times, method=method)
    print(f"{method.upper():12s} 95% CI: [{lower:.1f}, {upper:.1f}]")