Хабр, привет! Сегодня рассмотрим кейс, в котором классические статистические критерии не работают, и разберёмся, почему так происходит. Научимся строить свои собственные критерии по историческим данным. Обсудим плюсы и минусы такого подхода.

Меня зовут Коля, я работаю аналитиком данных в X5 Tech. Мы с Сашей продолжаем писать серию статей по А/Б тестированию. Это наша седьмая статья, предыдущие статьи можно найти в описании профиля.

Маленькие выборки из скошенных распределений

В этой статье будем рассматривать эксперименты на маленьких выборках из скошенных распределений. Скошенное распределение — это асимметричное распределение, у которого один из хвостов более тяжёлый. Пример скошенного распределения — логнормальное распределение. Формула и графики его функции плотности распределения представлены ниже.

Маленькая выборка — понятие относительное. Дадим неформальное определение, которое будем использовать в рамках этой статьи. Будем называть выборку маленькой, если распределение её среднего также имеет скошенное распределение.

Распределение среднего зависит от размера выборки. Если сгенерировать выборку размера 1 из скошенного распределения и посчитать среднее, то это среднее будет иметь такое же скошенное распределение. Если сгенерировать выборку огромного размера из скошенного распределения, то по ЦПТ распределение среднего будем стремиться к нормальному распределению. Продемонстрируем это на выборках размеров 1, 100 и 100 000. Данные будем генерировать из логнормального распределения с параметрами . На графиках ниже изображены оценки функций плотностей распределений среднего для разных размеров выборок.

Почему мы решили рассмотреть эксперименты на маленьких скошенных выборках? Оказывается, классические критерии, которые используются для проверки гипотезы о равенстве средних, работают некорректно на таких данных (про корректность А/Б тестов можно почитать в этой статье). Продемонстрируем это с помощью синтетических А/А тестов. Будем генерировать выборки размера 10 из логнормального распределения и проверять гипотезу о равенстве средних тестом Стьюдента. Повторим эту процедуру 10 000 раз и построим распределение p-value.

from collections import defaultdict

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [8, 6]

def plot_pvalue_distribution(dict_pvalues):
    """Рисует графики распределения p-value."""
    X = np.linspace(0, 1, 1000)
    for key, pvalues in dict_pvalues.items():
        Y = [np.mean(pvalues < x) for x in X]
        plt.plot(X, Y, label=key)
    plt.plot([0, 1], [0, 1], '--k', alpha=0.8)
    plt.title('Оценка распределения p-value', size=16)
    plt.xlabel('p-value', size=12)
    plt.legend(fontsize=12)
    plt.grid()
    plt.show()

sample_size = 10
sigma = 2
dict_pvalues = defaultdict(list)
for _ in range(10000):
    a, b = np.random.lognormal(sigma=sigma, size=(2, sample_size))
    _, pvalue = stats.ttest_ind(a, b)
    dict_pvalues['A/A'].append(pvalue)
plot_pvalue_distribution(dict_pvalues)

Распределение p-value оказалось неравномерным, критерий работает некорректно. При выводе теста Стьюдента используется предположение о нормальности распределения средних. Выше мы показали, что средние распределены не нормально. Вероятно, из-за этого критерий работает некорректно. Если увеличить размер выборок до 1 миллиона, то распределение средних будет близко к равномерному, и тест Стьюдента станет работать корректно.

Если в примере выше заменить критерий Стьюдента на критерий Манна-Уитни, то распределение p-value на синтетических А/А тестах станет равномерным. Однако, критерий Манна-Уитни в общем случае не проверяет гипотезу о равенстве средних. Это может приводить к противоречивым результатам. Например, для двух выборок с одинаковым средним критерий Манна-Уитни может отклонять нулевую гипотезу. Увеличение размера выборок не помогает решить эту проблему.

a = np.hstack([np.linspace(0, 1, 30), np.linspace(8, 9, 10)])
b = np.linspace(2, 3, 40)
pvalue = stats.mannwhitneyu(a, b).pvalue
print(f'разница средних = {a.mean() - b.mean()}')
print(f'pvalue = {pvalue:0.5f}')

разница средних = 0.0
pvalue = 0.00012

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

Как устроены статистические критерии

Мы предлагаем построить свой критерий для оценки подобных экспериментов. Сначала разберёмся, как вообще устроены статистические критерии. Вспомним ряд определений из математической статистики.

Выборка — набор случайных величин. В теории вероятностей и статистике отдельно выделяют случай независимых одинаково распределённых случайных величин (н.о.р.с.в.). Каждая из них распределена так же, как и остальные, и все они независимы в совокупности. Выборка из распределения F обозначается так: или .

Статистика — любая измеримая функция от выборки. Примеры статистик: минимум, среднее, квантиль. Статистику от выборки часто обозначают как .  

Статистическая гипотеза — любое предположение о распределении и свойствах случайной величины. Примеры статистических гипотез: выборка из нормального распределения, средние выборок равны. 

Уровень значимости — вероятность отклонить нулевую гипотезу при условии её истинности. Обозначается как .

Статистическая мощность — вероятность отклонения нулевой гипотезы в случае, когда альтернативная гипотеза верна.

Статистический критерий — математическое правило, позволяющее по реализациям выборок отвергнуть или не отвергнуть нулевую гипотезу с заданным уровнем значимости.

Критическая область — множество значений статистики критерия, при которых нулевая гипотеза отвергается.

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

• Если верна, то статистика имеет некоторое непрерывное распределение , в пределе не зависящее от ;

• Если неверна, то при .

Зная распределение статистики при верности нулевой гипотезы , можно определить область , в которую статистика попадает с вероятность . Конкретные границы области выбирают так, чтобы мощность критерия была максимальной. Если значение статистики попадает в область W, то нулевая гипотеза не отклоняется, иначе значение статистики попадает в критическую область, и нулевая гипотеза отклоняется.

Ниже на картинке приведён пример распределения G и области W для теста Стьюдента с .

Чтобы применять критерий, нужно знать распределение . Для этого математики доказывают теоремы, что для определённой гипотезы при некоторых предположениях распределение имеет такую-то функцию распределения. Если условия теоремы не выполняются, то критерий может работать некорректно.

Строим свой критерий

Чтобы построить критерий для проверки гипотезы о равенстве средних, нужно определить статистику и узнать её распределение при верности нулевой гипотезы. В качестве статистики логично использовать разность средних. Оценить распределение статистики можно по историческим данным с помощью синтетических А/А тестов. Для каждого синтетического А/А теста вычислим и запомним значение статистики, повторим это большое количество раз. По полученным значениям можно оценить функцию распределения. Функция распределения сама по себе для принятия решения нам не нужна, нам нужно p-value. Сравнив p-value с уровнем значимости, решаем, отклонять нулевую гипотезу или нет.

P-value — это вероятность получить заданное значение статистики или более экстремальное при верности нулевой гипотезы. Значение p-value вычисляется по формуле:

Чтобы оценить значение p-value двусторонней гипотезы, для реализации статистики t вычислим долю значений статистики из синтетических А/А тестов, которые оказались не больше t, и долю значений статистики из синтетических А/А тестов, которые оказались не меньше t. Возьмём меньшую из долей и умножим на два.

Реализуем описанный подход в коде. Допустим, данные имеют логнормальное распределение. Хотим проверить гипотезу о равенстве средних на группах размера 10. Есть 1 миллион исторических данных. Оценим распределение статистики по историческим данным с помощью 100 000 синтетических А/А тестов.

import seaborn as sns

def generate_statistic_values(history, sample_size, statistic_size):
    """Генерирует значения разности средних по историческим данным.

    history - исторические данные
    sample_size - размер выборки
    statistic_size - количество генерируемых данных
    
    return - массив со значениями разности средних
    """
    # для ускорения вычислений генерируем несколько групп за раз
    count_exp = len(history) // (2 * sample_size)
    statistic_list = []
    while len(statistic_list) < statistic_size:
        history_sample = np.random.choice(
            history,
            size=(count_exp, 2 * sample_size),
            replace=False
        )
        a_means = history_sample[:, :sample_size].mean(axis=1)
        b_means = history_sample[:, sample_size:].mean(axis=1)
        delta_means = b_means - a_means
        statistic_list += list(delta_means)
    statistic_array = np.array(statistic_list[:statistic_size])
    return statistic_array

alpha = 0.05
sample_size = 10
statistic_size = 100000
history_size = 1000000

history = np.random.lognormal(sigma=2, size=history_size)
statistic_values = generate_statistic_values(history, sample_size, statistic_size)

left_q, right_q = np.quantile(statistic_values, [alpha/2, 1-alpha/2])
kde = sns.kdeplot(statistic_values, label='оценка плотности \nраспределения G', gridsize=10000)
X, Y = kde.lines[0].get_data()
plt.plot([left_q, left_q], [-0.001, 0.006], 'k--', alpha=0.5)
plt.plot([right_q, right_q], [-0.001, 0.006], 'k--', alpha=0.5)
plt.plot([left_q, right_q], [0, 0], linewidth=4, label='область W')
plt.plot([min(X), max(X)], [0, 0], 'k')
for mask in [X < left_q, X > right_q]:
    plt.fill_between(X[mask], np.zeros(sum(mask)), Y[mask], color='k', alpha=0.1)
mask = (X > left_q) * (X < right_q)
plt.fill_between(X[mask], np.zeros(sum(mask)), Y[mask], color='yellow', alpha=0.1)
plt.plot([-30, -31], [0.0005, 0.009], 'k', alpha=0.5)
plt.plot([30, 31], [0.0005, 0.009], 'k', alpha=0.5)
for x in [-34, 29]:
    plt.text(x, 0.011, '$\\alpha / 2$', size=16)
plt.text(-4, 0.02, '$1-\\alpha$', size=16)
plt.legend(fontsize=12)
plt.xlim([-40, 40])
plt.show()

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

Проверим корректность полученного критерия. Проведём симуляцию 10 000 экспериментов с эффектом и без эффекта, построим распределения p-value.

def check_my_statistical_test(a, b, statistic_values):
    """Проверяет гипотезу о равенстве средних.

    a - значения контрольной группы
    b - значения экспериментальной группы
    statistic_values - значения разности средних

    return - pvalue
    """
    pe = b.mean() - a.mean()
    part_less = np.mean(statistic_values < pe)
    pvalue = 2 * min(part_less, 1 - part_less)
    return pvalue

effect = 30

dict_pvalues = defaultdict(list)
for _ in range(10000):
    a, b = np.random.lognormal(sigma=2, size=(2, sample_size))
    pvalue = check_my_statistical_test(a, b, statistic_values)
    dict_pvalues['A/A'].append(pvalue)
    b += effect
    pvalue = check_my_statistical_test(a, b, statistic_values)
    dict_pvalues['A/B'].append(pvalue)
plot_pvalue_distribution(dict_pvalues)

Критерий работает корректно, p-value на А/А тестах распределено равномерно, а на А/Б тестах выпукло вверх. 

Описанный подход построения критерия является весьма универсальным. Его можно использовать не только для проверки гипотезы о равенстве средних, но и для проверки других гипотез. Нужно выбрать подходящую статистику критерия и оценить её распределение по историческим данным при верности нулевой гипотезы.

Основным минусом подхода является необходимость иметь достаточно большое количество актуальных данных. Чем больше особенностей у распределения, тем больше потребуется исторических данных для построения корректного критерия. Исторические данные должны быть актуальные, то есть у них должно быть то же распределение, что и сейчас. Данные десятилетней давности, возможно, не подойдут.

Оценка распределения через синтетические А/А тесты – вычислительно затратная процедура. Если для вашей задачи классические критерии работают корректно и быстро, то разумно применять их.

Итоги

Если обычные критерии работают некорректно, можно построить свой критерий. Этот критерий будет заточен под определённую статистику и природу данных. Плюсы: работает для любых гипотез и данных. Минусы: нужны актуальные исторические данные, вычислительно затратно.

Алгоритм построения и применения критерия:

1. Выберите статистику критерия, подходящую для проверяемой гипотезы.

2. Сгенерируйте множество значений статистики с помощью синтетические А/А тестов на исторических данных.

3. Вычислите значение статистики на данных эксперимента.

4. Оцените значение p-value для реализации статистики из п. 3 по множеству значений из п. 2.