В этой статье мы рассмотрим, как оценить стандартное отклонение, построить доверительные интервалы и проверить гипотезы с помощью начальной загрузки. Мы узнаем, когда необходима начальная загрузка и каковы ее ограничения.
Метрики и точность их оценки
Представим, что мы работаем аналитиками в службе доставки интернет-магазина. Нам поставили задачу оценить, насколько быстро мы выполняем заказы. У нас есть данные о времени выполнения каждого заказа, и теперь нам нужно выбрать метрику и оценить ее значение. В этой статье предположим, что мы работаем с независимыми одинаково распределенными случайными величинами. Случай зависимых случайных величин будет рассмотрен в последующих статьях.
Самая простая метрика — среднее время выполнения заказа. Чтобы оценить среднее время выполнения заказа, мы можем взять все заказы за определенный период времени, например, за последний месяц, и рассчитать среднее время их выполнения.
Предположим, мы получили оценку среднего времени доставки, равного 90 минутам. Насколько надежна эта оценка? Понятно, что это, скорее всего, не истинное значение среднего времени доставки. Если мы подождем еще месяц и повторим расчет, мы получим чуть большее или чуть меньшее значение. Важно оценить стандартное отклонение полученной оценки, чтобы понять ее точность, поскольку 90±1 минута и 90±30 минут представляют совершенно разные ситуации.
Для среднего значения оценка стандартного отклонения рассчитывается по формуле:
где n — размер выборки, Xi — случайные переменные времени доставки, а Xavg — среднее значение выборки времени доставки.
Рассмотрим пример оценки среднего и стандартного отклонения. Допустим, у нас есть информация о 1000 доставок. Распределение времени доставки в реальной жизни может быть произвольным. В этом примере мы создадим время доставки из нормального распределения со средним значением 90 и стандартным отклонением 20. Мы создадим выборку, оценим среднее значение и оценим стандартное отклонение среднего значения.
import numpy as np
n = 1000
values = np.random.normal(90, 20, n)
mean = values.mean()
std = values.std() / np.sqrt(n)
print(f'Avg time to arrival estimation: {mean:0.2f}')
print(f'Std for avg time to arrival estimation: {std:0.2f}')
Оказалось, что в нашем примере достаточно 1000 наблюдений, чтобы стандартное отклонение средней оценки было меньше минуты.
квантили
Мы оценили среднее время выполнения заказа, что хорошо, но это всего лишь «средний показатель». Кто-то получает заказы быстрее, кто-то медленнее. Мы хотим, чтобы подавляющее большинство клиентов получали свои заказы достаточно быстро. Чтобы оценить, сколько времени требуется для доставки большинства заказов, мы можем использовать, например, 90-й процентиль. Каков физический смысл процентиля? Если 90-й процентиль равен 2 часам, это означает, что 90% заказов доставляются за 2 часа или меньше.
Мы можем легко оценить 90-й процентиль, используя данные, но как мы можем оценить стандартное отклонение этой оценки? Не существует простой универсальной теоретической формулы для оценки стандартного отклонения процентиля.
Было бы здорово, если бы у нас было 100 параллельных вселенных. В каждой вселенной мы собирали данные, вычисляли 100 процентилей и оценивали стандартное отклонение на основе полученных значений. Однако у нас нет 100 параллельных вселенных.
Кто-то может предложить разделить наши данные 1000 наблюдений на 10 частей по 100 значений в каждой. Затем вычислите значение квантиля для каждой части и оцените стандартное отклонение на основе этих 10 значений. Однако этот подход приведет к неверным результатам, поскольку стандартное отклонение оценки зависит от количества наблюдений, использованных для оценки значения квантиля. Чем больше у нас данных, тем точнее оценка и меньше стандартное отклонение.
Что нам делать тогда? Оказывается, есть метод, позволяющий оценить стандартное отклонение любой статистики, в том числе и квантиля. Это называется бутстрап.
Начальная загрузка
Bootstrap — это метод оценки стандартных отклонений и нахождения доверительных интервалов для статистических функционалов.
Давайте разберемся, как работает бутстрап. Помните, мы хотим оценить стандартное отклонение произвольной статистики. В этой статье мы оценим стандартное отклонение оценки 90-го процентиля.
Если бы мы могли получить данные из базового распределения, мы могли бы сгенерировать 100 выборок из этого распределения, рассчитать 100 квантилей на основе этих выборок и оценить стандартное отклонение. Однако мы не знаем истинного распределения, но можем оценить его по имеющимся данным.
В качестве оценки функции распределения будем использовать эмпирическую функцию распределения (ЭФР). EDF является несмещенной оценкой и сходится к истинной функции распределения по мере увеличения размера выборки.
Давайте посмотрим, как визуально выглядит эмпирическая функция распределения (EDF) для данных разных размеров из стандартного нормального распределения.
import matplotlib.pyplot as plt
from scipy import stats
def plot_ecdf(values, label, xlim):
X_ = sorted(set(values))
Y_ = [np.mean(values <= x) for x in X_]
X = [xlim[0]] + sum([[v, v] for v in X_], []) + [xlim[1]]
Y = [0, 0] + sum([[v, v] for v in Y_], [])
plt.plot(X, Y, label=label)
# data generation
for size in [20, 200]:
values = np.random.normal(size=size)
plot_ecdf(values, f'ЭФР, size={size}', [-3, 3])
# EDF
X = np.linspace(-3, 3, 1000)
Y = stats.norm.cdf(X)
plt.plot(X, Y, '--', color='k', label='ФР')
plt.legend()
plt.show()
График показывает, что по мере увеличения размера выборки эмпирическая функция распределения (EDF) лучше аппроксимирует истинную функцию распределения. Если мы увеличим размер выборки до нескольких тысяч, визуально отличить ФРЭ от истинной функции распределения станет сложно.
Поскольку мы знаем, что EDF обеспечивает хорошее приближение к истинной функции распределения, давайте сгенерируем из него данные! Как мы можем сделать это?
Генерация подвыборки размера n из EDF эквивалентна случайному выбору n элементов из исходной выборки с заменой. Это можно сделать всего одной строкой кода.
np.random.choice(values, size=n, replace=True)
Теперь мы можем оценить стандартное отклонение квантильной оценки. Для этого мы будем генерировать подвыборки из эмпирической функции распределения (EDF) 1000 раз. В каждой подвыборке мы рассчитаем значение статистики и оценим стандартное отклонение полученных значений.
n = 1000
B = 1000
values = np.random.normal(90, 20, n)
quantile = np.quantile(values, 0.9)
bootstrap_quantiles = []
for _ in range(B):
bootstrap_values = np.random.choice(values, n, True)
bootstrap_quantiles.append(np.quantile(bootstrap_values, 0.9))
std = np.std(bootstrap_quantiles)
print(f'90% quantile estimation: {quantile:0.2f}')
print(f'Std estimation for 90% quantile: {std:0.2f}')
Мы только что применили метод начальной загрузки для оценки стандартного отклонения 90-го процентиля.
Обратим внимание на два момента:
Чтобы обеспечить объективную оценку стандартного отклонения, необходимо создать подвыборки того же размера, что и исходная выборка.
Рекомендуется, чтобы количество итераций начальной загрузки находилось в диапазоне от 1000 до 10000. Обычно этого достаточно для получения достаточно точных результатов.
Доверительные интервалы
Мы научились оценивать стандартное отклонение статистики. Зная стандартное отклонение, можно интуитивно понять достоверность полученных результатов. Для получения более точных результатов можно построить доверительный интервал.
Доверительный интервал (ДИ) — это интервал, охватывающий оцениваемый параметр с заданной вероятностью.
Вернемся к примеру с доставкой из интернет-магазина и построим доверительный интервал для среднего времени доставки. В нашей статье мы будем использовать доверительный интервал 95%. Если у нас есть большой объем данных, независимо от распределения исходных данных (при условии, что дисперсия существует и не равна нулю), центральная предельная теорема утверждает, что среднее значение будет нормально распределено. Для нормально распределенной статистики доверительный интервал можно рассчитать по формуле:
где z(1−α/2) — квантиль стандартного нормального распределения, а α — уровень значимости. Для 95% доверительного интервала уровень значимости равен 0,05.
Чтобы построить доверительный интервал для квантиля с использованием метода начальной загрузки, необходимо выполнить три шага, аналогичных тому, как мы построили доверительный интервал для среднего значения:
- Оцените значение квантиля, используя исходные данные.
- Используйте начальную загрузку, чтобы оценить стандартное отклонение оценки квантиля.
- Рассчитайте доверительный интервал по формуле:
Такой доверительный интервал называется нормальным доверительным интервалом.
Нормальный доверительный интервал подходит, когда распределение статистики близко к нормальному распределению. Однако, если распределение статистики асимметрично, нормальный доверительный интервал может дать странные результаты. График ниже иллюстрирует ситуацию, когда граница доверительного интервала превышает минимальное значение распределения.
В случае асимметричных распределений мы можем использовать процентильный доверительный интервал. Чтобы построить процентиль CI, нам нужно отрезать площадь α/2 с каждой стороны распределения. Для 95% ДИ нам нужно отсечь 2,5%. На практике для расчета границ КИ нам необходимо оценить квантили qα/2 и q1−α/2 на основе статистика бутстрапа. Доверительный интервал будет иметь следующий вид:
Есть еще вариант — центральный доверительный интервал. Его границы определяются:
На первый взгляд этот доверительный интервал может показаться нелогичным, но он имеет строгое математическое обоснование.
Центральный доверительный интервал для логнормального распределения смещен в сторону нулевой плотности. Это можно интерпретировать как склонность к чрезмерной осторожности в случае, если у нас есть неполная информация о распределении, и на самом деле могут быть значения слева от нуля.
Мы узнали о трех методах построения доверительных интервалов с помощью начальной загрузки. Ниже приведен пример кода для построения этих интервалов.
import seaborn as sns
def get_normal_ci(bootstrap_stats, pe, alpha):
z = stats.norm.ppf(1 - alpha / 2)
se = np.std(bootstrap_stats)
left, right = pe - z * se, pe + z * se
return left, right
def get_percentile_ci(bootstrap_stats, pe, alpha):
left, right = np.quantile(bootstrap_stats, [alpha / 2, 1 - alpha / 2])
return left, right
def get_pivotal_ci(bootstrap_stats, pe, alpha):
left, right= 2 * pe - np.quantile(bootstrap_stats, [1 - alpha / 2, alpha / 2])
return left, right
n = 1000
B = 10000
alpha = 0.05
values = np.random.normal(90, 20, n)
quantile = np.quantile(values, 0.9)
bootstrap_quantiles = np.quantile(np.random.choice(values, (B, n), True), 0.9, axis=1)
normal_ci = get_normal_ci(bootstrap_quantiles, quantile, alpha)
percentile_ci = get_percentile_ci(bootstrap_quantiles, quantile, alpha)
pivotal_ci = get_pivotal_ci(bootstrap_quantiles, quantile, alpha)
sns.kdeplot(bootstrap_quantiles, label='kde stat')
plt.plot([quantile], [0], 'o', c='k', markersize=6, label='point estimation')
plt.plot([109, 120], [0, 0], 'k', linewidth=1)
d = 0.02
plt.plot(normal_ci, [-d, -d], label='normal CI')
plt.plot(percentile_ci, [-d*2, -d*2], label='percentile CI')
plt.plot(pivotal_ci, [-d*3, -d*3], label='central CI')
plt.legend()
plt.show()
Bootstrap и A/B-тестирование
Мы научились строить доверительные интервалы. Используя доверительный интервал, мы можем не только получить дополнительную информацию о значении метрики, но и проверить статистические гипотезы. Для проверки гипотезы о равенстве квантилей при уровне значимости 5 % достаточно построить 95 % доверительный интервал различия квантилей между группами. Если ноль находится за пределами доверительного интервала, то различия статистически значимы; в противном случае они не являются.
Опишем алгоритм проверки гипотез о равенстве двух произвольных метрик с помощью бутстрапа:
- Создайте пару подвыборок одинакового размера из исходных данных для контрольной и экспериментальной групп.
- Рассчитайте метрики (оценку метрик) для каждой группы.
- Вычислите разницу в метриках и сохраните полученное значение.
- Повторите шаги 1-3 от 1000 до 10000 раз.
- Постройте доверительный интервал с уровнем значимости α.
- Если 0 не принадлежит доверительному интервалу, то различия статистически значимы на уровне значимости α; в противном случае они не являются.
Давайте посмотрим, как этот алгоритм можно реализовать в коде. В качестве примера при формировании данных для экспериментальной группы уменьшим дисперсию, что приведет к уменьшению значения квантиля 90%.
n = 1000
B = 10000
alpha = 0.05
values_a = np.random.normal(90, 20, n)
values_b = np.random.normal(90, 15, n)
pe = np.quantile(values_b, 0.9) - np.quantile(values_a, 0.9)
bootstrap_values_a = np.random.choice(values_a, (B, n), True)
bootstrap_metrics_a = np.quantile(bootstrap_values_a, 0.9, axis=1)
bootstrap_values_b = np.random.choice(values_b, (B, n), True)
bootstrap_metrics_b = np.quantile(bootstrap_values_b, 0.9, axis=1)
bootstrap_stats = bootstrap_metrics_b - bootstrap_metrics_a
ci = get_percentile_ci(bootstrap_stats, pe, alpha)
has_effect = not (ci[0] < 0 < ci[1])
print(f'90% quantile value changed by: {pe:0.2f}')
print(f'{((1 - alpha) * 100)}% CI: ({ci[0]:0.2f}, {ci[1]:0.2f})')
print(f'Stat significant diff: {has_effect}')
Получается, что квантиль 90% значительно уменьшился.
В этом примере для построения доверительного интервала использовался процентильный доверительный интервал. На практике все три метода построения КИ часто дают близкие по точности результаты. Однако могут быть исключения, и результаты будут зависеть от характера данных и сравниваемых показателей. Чтобы определить, какой метод лучше всего работает в вашей ситуации, вам необходимо оценить ошибки первого и второго рода на основе исторических данных. Вы можете найти информацию о том, как это сделать, в нашей предыдущей статье о стратификации: как разделение выборки повышает чувствительность A/B-тестирования.
Ограничения Bootstrap
Bootstrap — очень универсальный метод, который позволяет численно исследовать свойства распределений произвольной статистики, а не только квантилей. Это делает его незаменимым для построения доверительных интервалов и проверки гипотез с нетривиальными метриками.
Если бутстрап так хорош, то почему он не используется во всех задачах? Основным недостатком бутстрапа является его вычислительная скорость. Применение бутстрапа — процедура, требующая значительных вычислительных ресурсов. Это становится заметно при работе с большими объемами данных и многократном применении бутстрапа. Расчеты могут занять часы, дни или даже недели.
Заметим, что во всех приведенных выше примерах предполагалось, что мы имеем дело с независимо и одинаково распределенными случайными величинами. В таких ситуациях случайные значения могут быть выбраны из исходной выборки для создания подвыборок. Однако если исходные данные зависимы и имеют сложный характер, это необходимо учитывать при формировании подвыборок. В противном случае оценка распределения статистики с использованием начальной загрузки может плохо аппроксимировать истинное распределение.
Важно не забывать о качестве исходных данных, используемых в бутстрапе. Если данные нерепрезентативны и плохо отражают истинное состояние, не следует ожидать надежных результатов.
Заключение
Мы обсудили статистический метод, который позволяет оценить стандартное отклонение и построить доверительные интервалы для нетривиальной статистики. Мы также обсудили, как применять бутстрап для проверки гипотез, где он незаменим, а также его ограничения.
