Бутстрап-методы для стандартного отклонения

Бутстрап: статистическая революция компьютерной эпохи

Бутстрап-ресамплинг — мощный статистический метод, который оценивает выборочное распределение любой статистики путём многократного ресамплинга из наблюдаемых данных. Предложенный Брэдли Эфроном в 1979 году, он произвёл революцию в статистическом выводе, позволяя анализировать сложные статистики без опоры на математические формулы или предположения о распределении.

Ключевая идея бутстрапа элегантно проста: ваша выборка — лучшая оценка генеральной совокупности. Ресамплируя из выборки (с возвращением), вы имитируете то, что происходило бы при многократном извлечении выборок из совокупности. Этот подход особенно ценен для стандартного отклонения, где традиционные формулы доверительных интервалов предполагают нормальность — допущение, которое часто нарушается на практике.

Бутстрап стал незаменимым в современной науке о данных, потому что работает с любой статистикой (медиана, корреляция, коэффициенты регрессии, веса нейросети) и не делает предположений о распределении данных.

Зачем бутстрап для стандартного отклонения?

Традиционные доверительные интервалы для стандартного отклонения предполагают, что данные взяты из нормального распределения. Когда это допущение нарушается (что случается часто), такие интервалы могут быть крайне неточными. Бутстрап предоставляет непараметрическую альтернативу.

Ключевые преимущества бутстрапа для стандартного отклонения:

• Никаких предположений о распределении: одинаково хорошо работает с нормальными, асимметричными и тяжелохвостыми данными
• Малые выборки: часто точнее параметрических методов при n < 30
• Сложные статистики: один и тот же подход применим к усечённому СО, MAD или любым мерам вариабельности
• Визуальное понимание: бутстрап-распределение показывает, что происходит, а не только итоговые числа

Процедура бутстрапа

Алгоритм бутстрапа удивительно прост. Из исходной выборки n наблюдений:

Сколько бутстрап-выборок? B = 1 000 часто достаточно для грубых оценок и проверки гипотез. Для доверительных интервалов B = 10 000 обеспечивает стабильные процентили. Для публикационного качества BCa-интервалов рекомендуется B = 15 000+.

Методы бутстрап-доверительных интервалов

Существует несколько методов построения доверительных интервалов по бутстрап-выборкам, каждый со своими компромиссами:

1. Процентильный метод (простейший)

Наиболее интуитивный подход: берутся процентили бутстрап-распределения напрямую.

Для 10 000 бутстрап-выборок это 250-е и 9750-е упорядоченные значения. Прост, но может быть смещённым при асимметричном бутстрап-распределении.

2. Базовый (стержневой) бутстрап

Использует связь между выборочной статистикой и бутстрап-статистиками:

Где θ̂ — исходное выборочное СО. Этот метод «отражает» процентильный интервал относительно выборочной оценки.

3. BCa (с коррекцией смещения и ускорением)

Золотой стандарт точности. BCa корректирует как смещение в бутстрап-распределении, так и ускорение (как стандартная ошибка меняется с параметром). Сложнее в расчёте, но обеспечивает интервалы второго порядка точности.

Решённый пример: ненормальные данные

Рассмотрим 15 измерений времени отклика (в мс): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312. Данные правосторонне-асимметричны (есть очень медленные ответы).

Бутстрап-ДИ асимметричен (шире с правой стороны), отражая правостороннюю асимметрию данных. ДИ по хи-квадрат не улавливает эту асимметрию.

Реализация на Python

Полная реализация бутстрапа с несколькими методами ДИ:

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