Metody bootstrapowe dla odchylenia standardowego

Bootstrap: rewolucja statystyczna ery komputerowej

Resampling bootstrapowy to potężna technika statystyczna, która estymuje rozkład próbkowania dowolnej statystyki poprzez wielokrotne ponowne próbkowanie z obserwowanych danych. Wprowadzony przez Bradleya Efrona w 1979 roku, zrewolucjonizował wnioskowanie statystyczne, umożliwiając analizę złożonych statystyk bez polegania na wzorach matematycznych czy założeniach o rozkładzie.

Kluczowa idea bootstrapu jest elegancko prosta: Twoja próbka jest najlepszym oszacowaniem populacji. Poprzez ponowne próbkowanie z Twojej próbki (ze zwracaniem) symulujesz, co by się wydarzyło, gdybyś mógł wielokrotnie próbkować z populacji. Podejście to jest szczególnie cenne dla odchylenia standardowego, gdzie tradycyjne wzory na przedziały ufności zakładają normalność — założenie, które często zawodzi w praktyce.

Bootstrap stał się niezbędny we współczesnej nauce o danych, ponieważ działa z dowolną statystyką (mediana, korelacja, współczynniki regresji, wagi sieci neuronowej) i nie zakłada niczego o rozkładzie danych.

Dlaczego bootstrap dla odchylenia standardowego?

Tradycyjne przedziały ufności dla odchylenia standardowego zakładają, że dane pochodzą z rozkładu normalnego. Gdy to założenie zawodzi (co jest częste), przedziały te mogą być dramatycznie niedokładne. Bootstrap zapewnia alternatywę niezależną od rozkładu.

Kluczowe zalety bootstrapu dla odchylenia standardowego:

• Brak założeń o rozkładzie: Działa równie dobrze z danymi normalnymi, skośnymi czy gruboogonowymi
• Wydajność przy małych próbkach: Często dokładniejszy niż metody parametryczne przy n < 30
• Obsługa złożonych statystyk: To samo podejście działa dla odchylenia uciętego, MAD czy niestandardowych miar zmienności
• Wizualny wgląd: Rozkład bootstrapowy pokazuje, co się dzieje, nie tylko końcowe liczby

Procedura bootstrapowa

Algorytm bootstrapowy jest niezwykle prosty. Z oryginalnej próbki n obserwacji:

Ile próbek bootstrapowych? B = 1 000 zwykle wystarcza do przybliżonych oszacowań i testów hipotez. Dla przedziałów ufności B = 10 000 zapewnia stabilne percentyle. Dla przedziałów BCa klasy publikacyjnej zaleca się B = 15 000+.

Bootstrapowe metody przedziałów ufności

Istnieje kilka metod konstruowania przedziałów ufności z próbek bootstrapowych, każda z kompromisami:

1. Metoda percentylowa (najprostsza)

Najbardziej intuicyjne podejście: weź percentyle rozkładu bootstrapowego bezpośrednio.

Dla 10 000 próbek bootstrapowych to 250. i 9 750. wartość uporządkowana. Proste, ale może być obciążone, gdy rozkład bootstrapowy jest skośny.

2. Podstawowy (pivotalny) bootstrap

Wykorzystuje zależność między statystyką próbkową a statystykami bootstrapowymi:

Gdzie θ̂ to oryginalne odchylenie próbkowe. Ta metoda “odbija” przedział percentylowy wokół oszacowania próbkowego.

3. BCa (z korektą obciążenia i przyspieszeniem)

Złoty standard pod względem dokładności. BCa koryguje zarówno obciążenie w rozkładzie bootstrapowym, jak i przyspieszenie (jak błąd standardowy zmienia się wraz z wartością parametru). Bardziej złożona obliczeniowo, ale zapewnia przedziały o dokładności drugiego rzędu.

Rozwiązany przykład: dane nienormalne

Rozważ 15 pomiarów czasu reakcji (w ms): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312. Dane są prawostronnie skośne (niektóre bardzo wolne odpowiedzi).

Bootstrapowy CI jest asymetryczny (szerszy w stronę wyższych wartości), odzwierciedlając prawostronną skośność danych. CI oparty na chi-kwadrat nie uchwytuje tej asymetrii.

Implementacja w Pythonie

Kompletna implementacja bootstrapu z wieloma metodami CI:

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}]")