Вам надоело каждый раз разбираться какую гипотезу, а главное с какими ограничениями к имеющимся данным проверяет бесчисленное множество статистических тестов?
Тогда бутстрап — это ваш выбор. Он не требует никаких параметрических предположений о данных или какой-либо нетривиальной математики и, вместе с тем, может быть применен к широкому спектру статистических оценок.
Вступление
Бутстрап позволяет из выборок, полученных из базы данных или в результате A/B-теста, путем повторного отбора наблюдений строить эмпирическое распределение любой выборочной статистики (метрики) без предварительных ограничений или требований к данным. С его помощью можно:
Построить доверительный интервал, например, для 60-го перцентиля, суммы или даже дисперсии;
Посчитать результаты эксперимента для медианы;
Найти p-value для приемочной ratio-метрики с зависимыми наблюдениями, такой как CTR или средний чек;
Провести анализ мощности разных статкритериев в A/B-тесте;
Постараться найти в какой части распределения произошел эффект от экспериментального воздействия.
❗️Основная идея❗️
Бутстрап старается притворяться генеральной совокупностью. Идейно можно к этому относиться так: у вас есть возможность провести сколько угодно «экспериментов» для проверки одной гипотезы, где у вас есть доступ ко всей «генеральной совокупности», которой является наша исходная выборка. Данная иллюзия поддерживается именно семплированием наблюдений с повторением.
Однако бутстрап не генерит новой информации и соответственно не повышает, репрезентативность исходных данных.
Кроме того, бутстрапированные выборки семплируются в объеме исходной. Это нужно, чтобы получить достаточно точную оценку вариативности (дисперсии) интересующей статистики, которая как раз и зависит от размера исходной выборки, как standard error для разницы средних в t-тесте.
Обобщенная реализация функции бутстрапа выглядит достаточно просто.
Взять исходную выборку;
Провести «эксперимент», отбирая в бут-выборку наблюдения с повторением;
На получившейся выборке посчитать интересующую статистику и положить ее в массив;
Повторить пункты 2 и 3 много-много раз и получить эмпирическое распределение статистики;
Получить из этого распределения нужную информацию;
Визуализировать для информативности.
import numpy as np
import pandas as pd
from scipy import stats
import seaborn as sns
import matplotlib.pyplot as plt
def bootstrap(sample, n_trials=5_000, statistic=np.median):
rng = np.random.default_rng()
stat_distrib = []
# get booted samples and count stat
for _ in range(n_trials):
boot_sample = rng.choice(sample, size=len(sample), replace=True)
stat_distrib.append(statistic(boot_sample))
result = do_some_math(stat_distrib)
do_some_viz(result, stat_distrib)
return result
def do_some_math(data):
pass
def do_some_viz(data):
passПрактика
Разберем указанные возможности бутстрапа на примерах. За основу возьмем набор непрерывных значений из экспоненциального распределения, например, трат пользователей. Данное распределение может выглядеть примерно так.
1. Построение доверительных интервалов
Давайте для нашей выборки значение 60-го перцентиля зададим интервальной оценкой с уровнем доверия 95%.
CONF_LVL = 0.95
def bootstrap_ci(sample, n_trials=5_000, statistic=np.median):
rng = np.random.default_rng()
stat_distrib = []
# get booted samples and count stat
for _ in range(n_trials):
boot_sample = rng.choice(sample, len(sample), replace=True)
stat_distrib.append(statistic(boot_sample))
result = do_some_math(stat_distrib)
do_some_viz(result, stat_distrib)
return result
def do_some_math(data):
# confidence interval counting
left_q = (1 - CONF_LVL) / 2
right_q = 1 - left_q
ci = np.quantile(data, [left_q, right_q])
return ci
def do_some_viz(res, data):
hist = plt.hist(data, bins=32, color='lightsalmon')
ymax = hist[0][np.argmax(hist[0])]
plt.vlines(np.mean(data), ymin=0, ymax=ymax, colors='black',
label=f'Statistic mean: {np.mean(data).round(3)}')
plt.vlines(res, ymin=0, ymax=ymax//2.5, linestyle='--', colors='black',
label=f'CI: {res[0].round(3)} and {res[1].round(3)}')
plt.xlabel('statistic value', fontsize=14)
plt.ylabel('frequency', fontsize=14)
plt.legend(loc=0)
return None
def quant_60(data):
return np.quantile(data, 0.6)
ci = bootstrap_ci(spendings, statistic=quant_60)На выходе получится приблизительно вот такой график эмпирического распределения выбранной статистики. Среднее значение данного распределения будет очень близко к оценке 60-го перцентиля в нашей исходной выборке.
На практике такое можно использовать для построения подневных графиков интересующих метрик с их границами ДИ или для презентации важной бизнесовой метрики, которая получилась в тестовой группе эксперимента, обычно это касается денег.
2. Результаты экспериментов по произвольной статистике
Иногда возникает продуктовая потребность оценить результаты эксперимента не для среднего значения, а для медианы. Реализовать это можно следующим образом.
CONF_LVL = 0.95
def bootstrap_ab(s1, s2, n_trials=5_000, statistic=np.median):
rng = np.random.default_rng()
stat_distrib = []
for _ in range(n_trials):
boot_s1 = rng.choice(s1, len(s1))
boot_s2 = rng.choice(s2, len(s2))
stat_distrib.append(statistic(boot_s1) - statistic(boot_s2))
result = do_some_math(stat_distrib)
do_some_viz(result[1], stat_distrib)
return result
def do_some_math(data):
# confidence interval counting
left_q = (1 - CONF_LVL) / 2
right_q = 1 - left_q
ci = np.quantile(data, [left_q, right_q])
# p_value
quant = stats.norm.cdf(x=0, loc=np.mean(data), scale=np.std(data, ddof=1))
p_value = quant * 2 if 0 < np.mean(data) else (1 - quant) * 2
return p_value, ci
p_value, ci = bootstrap_ab(spendings_t, spendings_c)Тут уже необходимо делать бут-выборки из двух «генеральных совокупностей», теста и контроля. Получив эмпирическое распределение разнозностей статистик, к нему можно относиться как в реализации «ЦПТ».
Чтобы найти p-value для двухсторонней гипотезы, необходимо узнать долю случаев, когда разница была равна 0 и случаев еще более выраженных. Для этого просто находим, в каком квантиле нашего эмпирического распределения расположен 0, и считаем хвосты.
Также в качестве критерия можно использовать ДИ. Если 0 за пределами границ ДИ, то отвергаем нулевую гипотезу H0.
Для проверки гипотезы о разницы средних значений существует разработанный статаппарат, а именно всеми любимый t-тест. В нем, согласно ЦПТ, разница средних имеет нормальное распределение, и для этого распределения известна мера шума (standard error) и оценки доверительных интервалов.
Так вот, бутстрап будет давать приблизительно такие же результаты при подсчете средних значений в эксперименте что и t-тест. Точность этих результатов будет расти вместе с ростом объема исходных выборок и количества итераций. Можно сказать, что бутстрап на средних значениях является вычислительно затратной аппроксимацией t-теста. Они также будут иметь приблизительно одинаковую мощность, но об этом ниже.
На примере со средними можно нагляднее продемонстрировать, почему бут-выборки формируются в объеме исходных данных. Если увеличить размер семплируемых бут-выборок в 25 раз, то эмпирическое распределение разниц среднего сузится, а standard deviation станет примерно в 5 раз меньше (если размеры групп приблизительно одинаковые) посчитанного значения standard error для t-теста.
3. Работа с ratio-метриками
Если мы для каждого пользователя поделим его траты на число заказов, то мы получим биномиальное распределение. Однако среднее значение по такому распределению будет средний чек на пользователя, т.е. среднее средних значений.
Обычно же средний чек в продукте определяется как сумма всех трат пользователей разделенная на сумму всех заказов. Такие показатели называются ratio-метриками, к ним еще относится CTR.
На картинке внизу представлено расхождение в значениях среднего чека на пользователя и среднего чека.
Суть проблемы ratio-метрик заключается в их оценки статистическими критериями, которые обычно требуют независимости входных величин. То есть, если построить распределение, где одно наблюдение это стоимость чека, то уже в нем всегда будут зависимые величины, так как несколько заказов мог сделать один пользователь, отсюда и зависимость.
Но и тут для бутстрапа нет никаких преград.
def bootstrap_ratio(s1_num, s1_denom, s2_num, s2_denom, n_trials=5_000):
rng = np.random.default_rng()
stat_distrib = []
for _ in range(n_trials):
users_idx1 = rng.choice(np.arange(0, len(s1_num)), size=len(s1_num))
boot1_num = s1_num[users_idx1]
boot1_denom = s1_denom[users_idx1]
ratio1 = boot1_num.sum() / boot1_denom.sum()
users_idx2 = rng.choice(np.arange(0, len(s2_num)), size=len(s2_num))
boot2_num = s2_num[users_idx2]
boot2_denom = s2_denom[users_idx2]
ratio2 = boot2_num.sum() / boot2_denom.sum()
stat_distrib.append(ratio1 - ratio2)
result = do_some_math(stat_distrib)
do_some_viz(result[1], stat_distrib)
return result
p_value, ci = bootstrap_ratio(spending_t, n_check_t, spending_c, n_check_c)В данной реализации надо семплировать с повторением случайных пользователей и отбирать их соответствующие сигналы для числителя (numerator) и знаменателя (denominator), а затем считать значение метрики. Построить эмпирическое распределение и получить из него ответы – это уже этап усвоенный.
Из картинки внизу видно, что среднее значение в распределении уже больше похоже на наблюдаемую разницу средних чеков.
4. Анализ мощности
Под мощностью теста обычно подразумевается вероятность принятия альтернативной гипотезы H1, когда она действительно верна. Мощность обычно зависит от размера имеющихся данных, а также модели эффекта, получаемого от экспериментального воздействия, и возможности статтеста этот эффект детектировать.
Поскольку у нас есть возможность провести сколько угодно много «экспериментов», можно взять и буквально посчитать долю событий, когда мы приняли H1 при заданном уровне alpha на разных статтестах. Возьмем t-test, mannwhitney test и bootstrap на средних.
ALPHA = 0.05
def power_compute(s1, s2, n_trials=200):
rng = np.random.default_rng()
# dict for stattests' p_values
test = {'t': [], 'mw': [], 'btsp': []}
for _ in range(n_trials):
boot_s1 = rng.choice(s1, len(s1))
boot_s2 = rng.choice(s2, len(s2))
test['t'].append(stats.ttest_ind(boot_s1, boot_s2)[1])
test['mw'].append(stats.mannwhitneyu(boot_s1, boot_s2)[1])
test['btsp'].append(bootstrap_ab(boot_s1, boot_s2,
n_trials=500, statistic=np.mean))
result = do_some_math(test)
do_some_viz(test)
return result
def do_some_math(test_dict):
power_dict = {}
for test in test_dict:
np_arr = np.array(test_dict[test])
power = len(np_arr[np_arr<=ALPHA]) / len(np_arr)
power_dict[f'{test}_power'] = power
return power_dict
def do_some_viz(test_dict):
df = pd.DataFrame(test_dict)
df = df.melt(var_name='test', value_name='p_val')
sns.ecdfplot(data=df, x='p_val', hue='test')
plt.vlines(x=ALPHA, ymin=0, ymax=1, colors='black')
plt.xlabel('Alpha')
plt.ylabel('Power')
return None
power_dict = power_compute(spendings_c, spendings_t)На выходе получится словарь со значениями мощностей для выбранных тестов. На картинке получится эмпирический кумулятивный (накопительный) график распределения (Empirical Cumulative Distribution Function). К нему можно относиться как к обычной гистограмме, в которой бины складываем слева направо, друг на друга. Однако ECDF более информативен, поскольку нам нужны не частоты, а суммарная доля случаев, когда p-value был ниже alpha. Доля интересующих нас случаев находится на пересечении вертикальной линии с линией соответствующего теста.
Из графиков видно, что мощность t-теста такая же как для бутстрапа на средних и находится на уровне 80%.
Анализ мощности с помощью бутстрапа может быть полезен для подведения результатов эксперимента.
Например, могут быть случаи, когда заранее остановили A/B-тест, потому что он прокрасился в положительную сторону, но после анализа мощности, выяснится, что корректно принять H1 можно на уровне 50%, что сравнимо c подбрасыванием монетки.
5. Перцентиль за перцентилем
Иногда в результатах A/B-тестов можно получать удивительные результаты, которые расходятся с нашими ожиданиями. Предположим, что был проведен эксперимент, и целевая метрика на мобилках статзначимо выросла, а на десктопе результат получился отрицательный, но серый. В качестве критерия использовался t-тест. Однако очень трудно пройти мимо такого падения и для десктопа можно узнать что непараметрический тест красит результат с мощностью 90%.
Поскольку тест Манна-Уитни работает с рангами, можно с помощью бутстрапа «пробежаться» по децилям тестовой и контрольной групп и найти, где произошло смещение рангов (позиций). В представленном ниже варианте это произошло в первых четырех децилях.
Итог
Понято, что бутстрап плохо масштабируется из-за своей вычислительной требовательности, особенно если в компании проведение A/B-тестов поставлено на поток, и за раз могут считаться десятки и сотни продуктовых метрик.
Однако для решения специфичных задач или при вычислении в экспериментах статзначимости самых неожиданных функций всегда можно попробовать «стукнуть» по своим данным грубой вычислительной палкой, чтобы они попробовали дать чуть больше информации. Если усвоить главную идею, которая стоит за реализацией бустрапа, все за этим следующее будет ограничиваться только фантазией применяющего его аналитика.
