Доверительные интервалы и бутстрэп

Оценка без интервала погрешности провоцирует слишком смелые выводы: модель с долей правильных ответов (accuracy) 85.9% лучше модели с 83.3% — если только тестовая выборка не мала. Заметка о байесовском выводе ответила байесовским доверительным интервалом (credible interval) — вероятностным утверждением о параметре, для которого нужно априорное распределение. Частотный (frequentist) ответ — доверительный интервал (confidence interval): процедура, которая на повторных выборках накрывает истинное значение в заданной доле случаев. В этой заметке мы объясняем, что означает эта доля, строим интервалы по центральной предельной теореме и вводим бутстрэп (bootstrap), который строит их по одним только данным, когда формулы нет.

Обозначения. $\theta$ — неизвестный параметр, $\hat\theta$ — его оценка, $\operatorname{se}$ — её стандартная ошибка. $z_{0.975} = 1.96$ — квантиль уровня 97.5% распределения $N(0, 1)$, а $t_{k, 0.975}$ — такая же квантиль $t$-распределения Стьюдента с $k$ степенями свободы.

Что означают 95%#

Доверительный интервал уровня 95% — это пара статистик $[L, U]$, функций от данных, таких, что

$$P_\theta\big(L \le \theta \le U\big) \ge 0.95 \quad \text{для любого } \theta .$$

Интервал с такой гарантией — при каждом $\theta$ и любом объёме выборки — называется точным, или консервативным, если покрытие может превышать 95%. Большинство интервалов на практике, в том числе интервалы нормального приближения, Уилсона и бутстрэпа ниже, — лишь номинальные 95-процентные интервалы: их покрытие стремится к 95% с ростом $n$ при условиях, на которых держится приближение, но при данных $n$ и $\theta$ оно может быть и меньше, и больше. Покрытие, полученное симуляцией при одном значении параметра, проверяет интервал в этой точке, а не при всех $\theta$.

Вероятность берётся по данным: интервал случаен, а параметр фиксирован. Если повторять эксперимент много раз, доля интервалов, содержащих истинное значение, приближается к их покрытию — здесь ровно 95% для t-интервала по нормальным данным:

Рисунок

Сто доверительных интервалов для одного и того же среднего

Когда данные получены, конкретный интервал либо содержит $\theta$, либо нет; утверждение «$\theta$ лежит в $[4.1, 6.3]$ с вероятностью 0.95» — это байесовское прочтение, и права на него доверительный интервал не даёт. На практике при больших выборках и слабых априорных распределениях интервалы двух видов почти совпадают, как показала заметка о байесовском выводе; различаются они, когда данных мало или априорное знание сильно.

Интервалы из нормального приближения#

Если оценка приближённо нормальна, $\hat\theta \approx N(\theta, \operatorname{se}^2)$, то интервал

$$\hat\theta \pm 1.96\, \widehat{\operatorname{se}}$$

имеет покрытие (coverage), близкое к 95%. Для среднего стандартная ошибка равна $\sigma/\sqrt n$, и $\sigma$ приходится оценивать выборочным стандартным отклонением $s$. При малом числе наблюдений эта оценка сама зашумлена, и интервал нужно расширить: если данные нормальны, величина $(\bar x - \mu)/(s/\sqrt n)$ имеет $t$-распределение Стьюдента с $n - 1$ степенями свободы, и интервал

$$\bar x \pm t_{n-1,\, 0.975}\, \frac{s}{\sqrt n}$$

имеет покрытие ровно 95%. При $n = 10$ множитель равен $2.262$, а не $1.96$; с множителем $1.96$ покрытие таких интервалов падает до 92%. С ростом $n$ разница исчезает: $t_{99, 0.975} = 1.984$.

Доли. Учебный интервал для доли, $\hat p \pm 1.96\sqrt{\hat p(1 - \hat p)/n}$ (интервал Вальда, Wald interval), ненадёжен, когда $n$ мало или $p$ близко к 0 или 1. При $n = 20$ и $p = 0.1$ его фактическое покрытие — 88%, а не 95%: всякий раз, когда не наблюдается ни одного успеха, а это происходит с вероятностью $0.9^{20} = 0.12$, интервал вырождается в одну точку 0. Интервал Уилсона (Wilson interval), который получается, если решить неравенство $\lvert \hat p - p \rvert \le 1.96 \sqrt{p(1-p)/n}$ относительно $p$, а не подставлять $\hat p$, в той же постановке имеет покрытие 96% и лучше подходит по умолчанию для долей кликов и долей правильных ответов на малых тестовых выборках. Но и он приближённый: при $n = 20$ его покрытие падает до 85% при $p = 0.008$.

Бутстрэп#

Для стандартной ошибки многих статистик нет удобной формулы: медианы, коэффициента корреляции, выигрыша одной модели над другой, а также метрик классификаторов, таких как площадь под ROC-кривой (AUC) и F1-мера из заметки о метриках классификации. Бутстрэп оценивает её симуляцией. Неизвестное распределение данных заменяется эмпирическим распределением, которое помещает массу $1/n$ в каждое наблюдение, а сэмплирование из генеральной совокупности — повторной выборкой (resampling) из данных с возвращением:

1. Извлечь из выборки $n$ значений с возвращением.
2. Вычислить статистику на этой повторной выборке.
3. Повторить $B$ раз (обычно от 1000 до 10 000).
4. Стандартное отклонение $B$ полученных значений оценивает стандартную ошибку; их перцентили уровней 2.5% и 97.5% дают 95-процентный перцентильный интервал (percentile interval).

Рисунок

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

Возьмём медиану 50 значений из экспоненциального распределения: у одной выборки медиана равна $0.549$ (истинная медиана — $\ln 2 = 0.693$). Бутстрэп-оценка стандартной ошибки по этой одной выборке равна $0.136$; истинная стандартная ошибка, найденная по 100 000 новых выборок, — $0.140$. Перцентильный интервал равен $[0.415, 0.999]$ и накрывает истинную медиану. Бутстрэп делает это, не зная, что данные экспоненциальны.

Когда бутстрэп не работает#

Бутстрэп работает, когда статистика гладко зависит от распределения данных, а выборка достаточно велика. В нескольких хорошо известных случаях он не работает или требует изменений:

• Экстремумы. Максимум повторной выборки никогда не превышает максимума исходной равномерной выборки и в точности совпадает с ним примерно в 63% повторных выборок ($1 - (1 - 1/n)^n$). Бутстрэп-распределение максимума совсем не похоже на истинное.
• Малые выборки. При 10 наблюдениях эмпирическое распределение плохо заменяет генеральную совокупность, и перцентильные интервалы дают заниженное покрытие; интервалы BCa (bias-corrected and accelerated, с поправкой на смещение и ускорение) помогают, но ненамного.
• Зависимые данные. Повторная выборка отдельных наблюдений разрушает зависимость. Временным рядам нужен блочный бутстрэп (block bootstrap), а данным, разбитым на группы (несколько примеров на пользователя или на документ), — повторная выборка целых групп.
• Тяжёлые хвосты. Для выборочного среднего данных с бесконечной дисперсией, например с показателем хвоста меньше 2, бутстрэп несостоятелен.

Сравнение двух моделей на одной тестовой выборке#

Два классификатора оценены на одних и тех же 1000 примерах, их доли правильных ответов — $0.833$ и $0.859$: реальна ли разница? Ошибки двух моделей сильно коррелированы, потому что трудные примеры трудны для обеих. Парный (paired) бутстрэп повторно выбирает примеры, сохраняет вместе результаты обеих моделей на каждом выбранном примере и вычисляет разницу на каждой повторной выборке. Его 95-процентный интервал для разницы, $[0.006, 0.046]$, не содержит нуля. Если считать две доли независимыми и выбирать их повторно по отдельности, получится $[-0.005, 0.057]$, и этот интервал ноль содержит: игнорируя парность, мы выбрасываем информацию о том, что модели ошибаются на одних и тех же примерах. Тот же принцип — сравнивать на одних и тех же единицах и повторно выбирать именно единицы — применим к A/B-тестам с повторными измерениями и к оценке моделей на наборах промптов или документов.

Что изменилось#

	Точечная оценка	Доверительный интервал
Сообщает	$\hat\theta$	$[L, U]$ с заявленным покрытием
Смысл 95%		при повторных выборках процедура накрывает $\theta$ не менее чем в 95% случаев
Среднее, малое $n$		$\bar x \pm t_{n-1} s/\sqrt n$
Доля	$\hat p$	интервал Уилсона, а не Вальда
Нет формулы		бутстрэп: повторная выборка, пересчёт, перцентили
Две модели, одна тестовая выборка	два числа	парный интервал для разницы

Ограничения#

• Покрытие — обещание о процедуре при её предположениях. Независимые наблюдения, верная модель, правильная стандартная ошибка; кластеризованные данные и подгонка анализа под данные (data snooping) незаметно его нарушают.
• Бутстрэп оценивает выборочную изменчивость, но не устраняет систематическую ошибку. Если выборка нерепрезентативна, повторная выборка из неё воспроизводит ту же ошибку; смещение оценки бутстрэп может оценить, но перцентильный интервал его не исправляет.
• Один интервал отвечает на один вопрос. Привести много интервалов и выделить те, что не содержат нуля, — это проблема множественных сравнений, предмет следующей заметки.

Что нужно запомнить#

• Точный доверительный интервал уровня 95% строится процедурой, которая накрывает истинное значение не менее чем в 95% повторных выборок при любом значении параметра; интервалы нормального приближения, Уилсона и бутстрэпа достигают 95% лишь приближённо. Каждый отдельный интервал либо накрывает значение, либо нет.
• $\hat\theta \pm 1.96\, \widehat{\operatorname{se}}$ для приближённо нормальных оценок; для среднего при малом $n$ — $\bar x \pm t_{n-1, 0.975}\, s/\sqrt n$.
• Интервал Вальда для доли даёт заниженное покрытие при малом $n$ или крайних $p$; лучше интервал Уилсона.
• Бутстрэп: повторная выборка с возвращением, пересчёт статистики, стандартное отклонение или перцентили полученных значений.
• Он не работает для экстремумов, очень малых выборок, зависимых данных (нужны блоки или группы) и бесконечной дисперсии.
• Две модели на одной тестовой выборке сравнивают парным интервалом.

Код#

Покрытие t-интервала и интервалов Вальда и Уилсона, бутстрэп медианы и максимума и парное сравнение:

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# coverage of 95% intervals for a mean, n = 10 from N(5, 2^2): the t quantile 2.262 against the normal 1.96
X = rng.normal(5, 2, (100_000, 10))
m, s = X.mean(axis=1), X.std(axis=1, ddof=1)
for q in (stats.t(9).ppf(0.975), 1.96):
    h = q * s / np.sqrt(10)
    print(round(float(q), 3), round(float(((m - h <= 5) & (5 <= m + h)).mean()), 3))
# 2.262 0.951 | 1.96 0.919

# a proportion with n = 20, p = 0.1: coverage of the Wald and the Wilson intervals
n, p, z = 20, 0.1, 1.96
ph = rng.binomial(n, p, 200_000) / n
w = z * np.sqrt(ph * (1 - ph) / n)
c = (ph + z * z / (2 * n)) / (1 + z * z / n)
h = z * np.sqrt(ph * (1 - ph) / n + z * z / (4 * n * n)) / (1 + z * z / n)
print(round(float(((ph - w <= p) & (p <= ph + w)).mean()), 3), round(float(((c - h <= p) & (p <= c + h)).mean()), 3))   # 0.877 0.957

# the bootstrap for the median of 50 exponential values: standard error and percentile interval
x = rng.exponential(1, 50)
boot = np.median(rng.choice(x, (2000, 50), replace=True), axis=1)
true = np.median(rng.exponential(1, (100_000, 50)), axis=1)
print(round(float(np.median(x)), 3), round(float(boot.std()), 3), np.percentile(boot, [2.5, 97.5]).round(3), round(float(true.std()), 3))
# 0.549 0.136 [0.415 0.999] 0.14

# the bootstrap fails for the maximum: the sample maximum reappears in about 63% of resamples
u = rng.random(50)
print(round(float((rng.choice(u, (2000, 50)).max(axis=1) == u.max()).mean()), 3), round(1 - (1 - 1 / 50) ** 50, 3))   # 0.625 0.636

# comparing two models on the same 1,000 test examples: paired against unpaired bootstrap
N = 1000
difficulty = rng.normal(size=N)
a = difficulty + 0.4 * rng.normal(size=N) > -1.04         # model A is correct on these examples
b = difficulty + 0.4 * rng.normal(size=N) > -1.18         # model B, correct mostly on the same ones
idx, idx2 = rng.integers(0, N, (5000, N)), rng.integers(0, N, (5000, N))
paired, unpaired = b[idx].mean(axis=1) - a[idx].mean(axis=1), b[idx2].mean(axis=1) - a[idx].mean(axis=1)
print(a.mean(), b.mean(), np.percentile(paired, [2.5, 97.5]).round(3), np.percentile(unpaired, [2.5, 97.5]).round(3))
# 0.833 0.859 [0.006 0.046] [-0.005  0.057]

Дальше#

Доверительный 95-процентный интервал для разницы, не содержащий нуля, отвергает гипотезу об отсутствии разницы на уровне 5% при предположениях, на которых построен интервал. Он не доказывает, что разница реальна или достаточно велика, чтобы иметь значение, и значит меньше, если был одним из многих просмотренных интервалов. Следующая заметка делает это решение явным: нулевые гипотезы, статистики критериев, p-значения и многочисленные ошибки в их прочтении, мощность — способность обнаружить эффект — и A/B-тест от планирования до ловушек подглядывания и множественных сравнений.

Литература#

• J. Neyman. Outline of a Theory of Statistical Estimation Based on the Classical Theory of Probability. Philosophical Transactions of the Royal Society A, 236, 1937.
• B. Efron. Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics, 7(1), 1979.
• B. Efron, R. J. Tibshirani. An Introduction to the Bootstrap. Chapman & Hall, 1993.
• L. D. Brown, T. T. Cai, A. DasGupta. Interval Estimation for a Binomial Proportion. Statistical Science, 16(2), 2001.
• L. Wasserman. All of Statistics. Springer, 2004. Глава 6 (доверительные интервалы), глава 8 (бутстрэп).