Решаем простую статистическую задачу пятью способами18.05.2024 16:30

Вот тут ув. @dimview на пальцах и Си объясняет за бутстрап решая несложную задачу. И в статистике существует 100500 разных тестов для (не)подтверждения нулевой гипотезы.

Давайте используем ряд самых распространеных и посмотрим на результаты. В конце сравним с бутстрапом. Изложение будет сопровождаться кратким выводом и объяснением основных тестов, их «ручной» реализацией и сравнением результата с готовыми тестами из пакета scipy.stats. В этом плане, мне кажется, повторение лишним не будет, т.к. позволит лучше понять и уяснить принцип и особенности тестов.

Сама задача звучит как: «И вот свежие результаты — в тестовой группе из 893 пришедших у нас что-то купили 34, а в контрольной группе из 923 пришедших что-то купили 28. Возникает вопрос — идти к начальству и говорить «в тестовой группе конверсия 3.81%, в контрольной группе 3.03%, налицо улучшение на 26%, где моя премия?» или продолжать сбор данных, потому что разница в 6 человек — ещё не статистика?»

Сравнение средних

Первый подход — рассмотреть продажи как 1, а остальные посещения как 0. И сравним конверсии (т.е. отношение положительных исходов к общему числу испытаний). И первое, приходит на ум — использовать двувыборочный t-тест.

Пусть число продаж в первом случае — k1 (=34), а всего пользователей n1 (=893). И k2, n2 соответственно для контрольной группы. Мы хотим посмотреть на распределение случайной величины k1/n1 — k2/n2 если бы эксперименты много раз проводились в параллельных мирах. Но при этом бы выполнялись условия: n1 и n2 были бы фиксированы, а k1 + k2 = const. Как раз в такой ситуации у нас возникает t-распределение.

Самое простое — сделать два массива и напрямую подставить в t-test из пакета scipy.stats и получить p-value.

import scipy.stats as ss
#Тестовая бырока
n1 = 893
k1 = 34
#Контрольная выборка
n2 = 923
k2 = 28

s1 = np.zeros(n1)
s1[:k1] = 1
s2 = np.zeros(n2)
s2[:k2] = 1
ss.ttest_ind(s1,s2)
>>TtestResult(statistic=0.9075373501154826, pvalue=0.36424324967099464, df=1814.0)
#s1 и s2 - это представление нашей выборки в виде 1 - продаж
#и 0 - их отсутствия. Это нам понадобится еще и на бутстрапе.
print(s1[:70])
>>[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

Напомню как устроен t-test. Мы смотрим на распределение разницы выборочных средних (двувыборочный t-критерий).

При сложении независимых случайных величин их дисперсии также складываются ( $D(\xi+\mu)=D\xi+D\mu$ ). Поделив нашу разность на стандартное отклонение получим t-распределение.

$t=\frac{ \overline{x_1}-\overline{x_2} } { \sqrt{ \frac{D_1}{n_1}+\frac{D_2}{n_2}} }$

Количество степеней свободы df = n1+n2–2.

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

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

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

Данные выражения получаются из элементарной арифметики если подставить их в формулы для несмещенной дисперсии

d1 = k1/(n1-1)*(1-k1/n1)
d2 = k2/(n2-1)*(1-k2/n2)
D = d1+d2

m = k1/n1 - k2/n2
std = np.sqrt( d1/n1 + d2/n2 )
t = m/std
print('t-value: ',t)
df = n1+n2-2
print('p-value: ',1 - ss.t.cdf(t, df) )
>>t-value:  0.905901136585428
>>p-value:  0.18255424508928209
# Наше p-value получилось почти в 2 раза меньше, чем из теста.
# Потому что мы посчитали односторонее значение. 
# Для двустороннего значения, строго говоря, нужно вычислить вот такую штуку:
print('two-side p-value: ',ss.t.cdf(-t, df) + 1 - ss.t.cdf(t, df) )
>>two-side p-value:  0.36510849017856417
#Или неспортивно домжножить на 2, т.к. распределение симметричное 
print('two-side p-value: ', (1 - ss.t.cdf(t, df))*2 )
two-side p-value:  0.36510849017856417
#Ну и ради интереса сравним с нормальным распределением
print('two-side p-value: ', (1 - ss.norm.cdf(m, scale=std) )*2 )
>>two-side p-value:  0.3649881826944741

Как видим, значения p-value для нормального и t-распределений отличаются фактически в 4 м знаке после запятой, т.к. количество степеней свободы оч. большое и можно не запариваться с t-распределением. Еще у теста из пакета scipy немного меньшее t-значение (не знаю, почему).

# Ради интереса доверительные интервалы 95% для разности средних
print(f'{ss.t.ppf(.025,df,loc=m,scale=std)*100:.3f}%; {ss.t.ppf(.975,df,loc=m,scale=std)*100:.3f}%')
>>-0.901%; 2.449%
# И для нормального
print(f'{ss.norm.ppf(.025,loc=m,scale=std)*100:.3f}%; {ss.norm.ppf(.975,loc=m,scale=std)*100:.3f}%')
>>-0.900%; 2.448%

Критерий

Ранее мы рассматривали задачу как случайную величину с некоторым средним. Теперь мы можем представить, что бросаем несимметричную монетку. Отличие ожидаемых частот оценивается расстоянием Пирсона $\chi^2 = \sum \frac{(O-E)^2}{E}$ , где O и E — observed and expected частоты. И это расстояние имеет $\chi^2$ распределение.

Например мы бросаем монетку 60 раз. Мы ожидаем (expect), что орел выпадет 30 раз, 60–30 — выпадет решка.

А наблюдаем (observe) 20 решек и 40 орлов. Должно быть E (орла)=Е (решка)=30. Получили О (решка)=20, О (орла)=40.

Отсюда

$\chi^2 = \left(\frac{20-30}{\sqrt{30}}\right)^2+\left(\frac{40-30}{\sqrt{30}}\right)^2$

Теперь вернемся к нашей задаче. У нас есть 2 измерения с частотами. Какова вероятность, что распределения частот между ними на самом деле одинаковы, а наблюдаемая разница случайна (принимаем H_0 — нулевую гипотезу)? Порешаем прошлый пример.

Rоличество степеней свободы: df = (m-1)(n-1)
В данном случае равно (2–1)(2–1) = 2

# Найдем общее среднее.
# Можно интерпретировать это как количество положительных
# исходов на общее число испытаний, если в контрольной
# и тестовой выборках чатота на самом деле одинакова.
total_mean = (k1+k2)/(n1+n2)

#Теперь найдем частоты, какими бы они должны были
#быть с учетом дисбаланса классов
k1_e = n1*total_mean
k2_e = n2*total_mean
print(f'k1 expected = {k1_e:.1f}; k2 expected = {k2_e:.1f}')
>>k1 expected = 30.5; k2 expected = 31.5

#Это для частот "промахов"
t1 = n1 - k1
t2 = n2 - k2

t1_e = n1*(t1+t2)/(n1+n2)
t2_e = n2*(t1+t2)/(n1+n2)

#Или можно просто использовать 1 - total_mean
#вместо (t1+t2)/(n1+n2)

print(f't1 expected = {t1_e:.1f}; t2 expected = {t2_e:.1f}')
>>t1 expected = 862.5; t2 expected = 891.5

O = np.array([k1,k2,t1,t2])
E = np.array([k1_e,k2_e,t1_e,t2_e])

chi2 = sum( (O-E)**2/E )
print('расстояние chi2: ',chi2)
>>расстояние chi2:  0.8241579182647449
print(f'p-value: {1 - ss.chi2(1).cdf(chi2):.3f}', )
>>p-value: 0.364

Заметим, что ничего тут на 2 не умножаем, т.к. распределение несимметричное и одностороннее. Ну и воспользуемся тестом $\chi^2$ из scipy

#Тут мы убрали поправку Йейтса
res = ss.chi2_contingency( [[k1,k2],[n1-k1,n2-k2]], correction=False)
print(f'Расстояние хи2: {res.statistic:.3f}')
>>Расстояние хи2: 0.824
print(f'p-value: {res.pvalue:.3f}')
>>p-value: 0.364

Как видим, результаты в точности совпадют с таковыми для t-теста. Что не очень удивительно, учитывая эквивалентность этих тестов для матрицы сопряженности размером 2×2.

ANOVA

Однофакторный дисперсионный анализ.

В первом случае мы сравнивали средние значения и получали t-распределение. Во втором — дисперсии и получали хи2 распределение. В случае с anova мы будем сравнить отношение дисперсий и получим, соответственно, F-распределение.

Anova основана на том, что есть несколько категориальных переменных, которые в среднем принимают различные значения. В нашем случае — это тестовая и контрольная группа посетителей на сайт. Есть изменчивость данных, которая складывается из изменчивости внутри групп так и изменчивости между группами (внутригрупповая и межгрупповая дисперсии). На цифрах это вот так выглядит.

Пусть есть группы x_i^j , где j=1…M — количество групп (в нашем случае 2), а i=1…N — количество элементов в группе (в нашем случае 893 и 923).

Сначала мы находим общегрупповое среднее

$\bar{\bar{x}} =\frac{1}{MN} \sum_{i,j} x_i^j$

Потом найдем общегрупповую сумму квадратов

$SST = \sum_{i,j} (x_i^j-\bar{\bar{x}})^2$

SST раскладывается на SST=SSW+SSB

Теперь рассчитаем внутрегрупповую сумму квадратов

$SSW = \sum_j \sum_i(x_i^j-\bar{x}^j)^2$

Т.е. для каждой группы берется свое среднее ${x}^j$ . $df = N\cdot M-M$ . Т.е. общее количество элементов минус количество групп.

Теперь межгруппову сумму квадратов

$SSB = \sum_j N_j( \bar{x}^j-\bar{\bar{x}} )^2$

Каждое слагаемое домножили на количество элементов в группе ( N_j ). Количество степеней свободы df = M-1 . Тут вроде понятно — знаем общегупповое среднее и пару средних — можем вычислить третье.

Теперь можно посчитать F значаение

$F = \tfrac{ \dfrac{SSB}{M-1} }{ \dfrac{SSW}{N\cdot M-M} }$

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

mean1 = k1/n1
mean2 = k2/n2
total_mean = (k1+k2)/(n1+n2)

ssw = k1*(1-mean1)**2 + (n1-k1)*mean1**2 + k2*(1-mean2)**2 + (n2-k2)*mean2**2
df_ssw = n1+n2-2

ssb = n1*(mean1-total_mean)**2 + n2*(mean2-total_mean)**2
df_ssb = 2 - 1

F = ssb/df_ssb / (ssw/df_ssw)
print('F-value: ',F )
#>>F-value:  0.8236240418546322
print('P-value: ',1 - ss.f.cdf(F,df_ssb,df_ssw) )
#>>P-value:  0.36424324967094446

Получили ровно такое же p-value.

В отличии от χ² anova может сравнивать непрерывные данные (т.е. ему не обязательно частоты вычислять). И может проводить множественные сравнения в отличии от t-теста.

Про χ² на Хабре есть обстоятельная статья.

Точный тест Фишера

Точный тест Фишера как бы «в лоб» подсчитывает вероятность получить такие или более выраженные отклонения наблюдаемых значений при условии верности нулевой гипотезы. В нашем случае это то, что подвыборки (k1, n1) и и (k2, n2) получились из одной выборки.
Т.е. у нас есть k1+k2 положительных исходов в N = n1+n2 испытаниях.

Давайте построим матрицу

$\begin{matrix}{} k_1 & k_2 & k_1 + k_2\\ \nu_1 & \nu_2 & \nu_1+\nu_2\\ n_1 & n_2 & N\end{matrix}$

Третья колонка и третий столбец фиксированный. Т.е. у нас фиксировано общее количество положительных/отрицательных исходов, количество экспериментов в кажом тесте ( $\nu_{1,2} = n_{1,2} - k_{1,2}$ ) и, соответственно, общее количество экспериментов (показов).

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

Первый шаг

Давайте построим матрицу

$\begin{matrix}{} k_1 & k_2 & k_1 + k_2\\ \nu_1 & \nu_2 & \nu_1+\nu_2\\ n_1 & n_2 & N\end{matrix}$

Третья колонка и третий столбец фиксированный. Т.е. у нас фиксировано общее количество положительных/отрицательных исходов, количество экспериментов в кажом тесте ( $\nu_{1,2} = n_{1,2} - k_{1,2}$ ) и, соответственно, общее количество экспериментов (показов).

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

Первый шаг

Рассмотрим наш «процесс» как ряд событий.

Событие С: К нам пришли N=n1+n2 пользователей и k1+k2 из них у нас что-то купили. Вероятность этого события посчитается как $C_N^{k_1+k_2}p^{k_1+k_2}(1-p)^{N}$ . Как видим — это выборка из Биномиального распределения с вероятностью продажи .

Событие А: В тестовой группе было k1 продаж из n1 пользователей

$\frac{n_1!}{k_1!(n_1-k_1)!}p^{k_1}(1-p)^{\nu_1} = C_{n_1}^{k_1}p^{k_1}(1-p)^{\nu_1}$

Событие B: То же самое, только для пары (k2, n2)

Теперь нужно найти вероятность получить событие AB при условии, что произошло C: P(AB|C) . Воспользуемся формулой Байеса и независимостью событий A и B.

$P(AB|C) = \frac{P(C|AB)P(AB)}{P(C)} = \frac{P(AB)}{P(C)} = \frac{P(A)P(B)}{P(C)}$

Тут применено соображение что P(C|AB)=1 т.к. AB являетс следствие C. Далее легко подставить наши вычисленные вероятности:

$P(AB|C) = \frac{ C_{n_1}^{k_1}p^{k_1}(1-p)^{\nu_1} C_{n_2}^{k_2}p^{k_2}(1-p)^{\nu_2} } {C_{n_1+n_2}^{k_1+k_2}p^{k_1+k_2}(1-p)^{\nu_1+\nu_2}} = \frac{ C_{n_1}^{k_1} C_{n_2}^{k_2} } {C_{n_1+n_2}^{k_1+k_2}}$$ Это не что иное как гипергеометрическое распределение с параметрами $$ HG(k_1;n_1+n_2,k_1+k_2,n_1) = \frac{ C_{k_1+k_2}^{k_1} C_{n_1+n_2-k_1-k_2}^{n_1-k_1} } { С_{n_1+n_2}^{n_1} }$

Убедиться в этом можно арифметически.

Мы меняем параметр k_1 на k_1' т.е. в первой группе было не 34 продажи, а 35 или 36 например. При этом должно быть соблюдено условие, где k_1 + k_2 = k_1' + k_2' . Т.е. если в первой группе было 35 продаж, то во второй — 27 вместо 28.

Второй шаг

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

#Для начала убедимся, что у нас действительно гипергеометрическое распределене
print( math.comb(n1,k1)*math.comb(n2,k2)/math.comb(n1+n2,k1+k2) )
#>>0.06830756250854167
print( ss.hypergeom(n1+n2,k1+k2,n1).pmf(k1) )
#>>0.06830756250854171
#Использование такой формы избавляет нас от того, чтобы прибавлять к k1 в первом comb и вычитать из k2 во втором.

Вероятность получить такое же или большее отклонение опишется как можно найти:

p_all = ss.hypergeom(n1+n2,k1+k2,n1).pmf(np.arange(0,63))
p_less = p_all[p_all <= ss.hypergeom(n1+n2,k1+k2,n1).pmf(34)]
print(f'''Двухсторонний тест: {p_less.sum():.4f} ''')
#>>Двухсторонний тест: 0.3694 
print(f'''Чрез функцию scipy: {ss.fisher_exact(np.array([[34,859],[28,895]])
                              ,alternative='two-sided')[1]:.4f} ''')
#>>Чрез функцию scipy: 0.3694 

# Ради интереса выведем аргументы k1, где вероятность такая же или меньше
np.array( [np.where(p_all==i)[0][0] for i in p_less] )
>>array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 34, 35, 36, 37, 38, 39, 40,
       41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57,
       58, 59, 60, 61, 62], dtype=int64)

Значение двухстороннего p-value почти точно совпадает с значениями двухсторонних тестов выше (хотя в тысячных различия есть). ANOVA c t-тестом и Хи2 фактически могут мерить только двухсторонние p-value (т.к. они симметричные/квадратичные).

Точный Фишер может и одностороннее значение померять, которое отличается на 0.035 сотых, что довольно не мало.

При «точном Фишере» не нужно особо забоититься о минимальных частотах, гауссовости выборочных средних. В общем там где считается быстро, данный тест выглядит привлекательно и легко интерпретируется.

Bootstrap для народа

Ну и давайте перейдем к, так сказать, численным методам. Вот мы получили две выборки. Давайте будем считать, что это и есть генеральная совокупность. Т.е. Данные у нас бинарные, можно представить их как 2 мешка с черными и белыми шарами, из которых мы вынимаем шары с возвращением. Соответсственно мы можем вынуть сколько угодно шаров. Это наша имитация генеральной совокупности.

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

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

import plotly.express as px
from scipy.stats import bootstrap

Мы считаем, что наши продажи имеют номера от 0 до k1–1. И мы генерируем слуайные числа в диапазоне от 0 до n1–1. Количество числе < k1 будет количеством положительных исходов (продаж) в данной выборке в тестовой группе. То же самое в контрольной.

Мы делаем M таких выборок с возвращением в тестовой и контрольной группе. Каждый раз считаем разницу конверсий в группах, т.е. k_1^m/n_1 - k_2^m/n_2 , где m=0...M-1 . Это такая «логическая имитация» биномиального процесса

#Создали дефолтный генератор
rng = np.random.default_rng()
M = 9999
dlt = np.zeros(M)
m = k1/n1 - k2/n2
for i in range(M):    
    smpl1 = rng.integers(n1,size=n1)
    k11 = (smpl1= m) | (dlt - m <= -m) ).sum()/M )
#p-value:  0.36503650365036505

Можно и напрямую нагенерировать выборок из биномиального распределения

sales1 = rng.binomial(n1,k1/n1,M)
sales2 = rng.binomial(n2,k2/n2,M)
dlt2 = sales1/n1 - sales2/n2
print('p-value: ', (  (dlt2 - m >= m) | (dlt2 - m <= -m) ).sum()/M )
#>>p-value:  0.35733573357335735

Или с помощью функции scipy

def difmean(s1,s2):
    return s1.sum()/s1.shape[0] - s2.sum()/s2.shape[0]
res = bootstrap((s1,s2), difmean, confidence_level=0.95, n_resamples=M, )
dlt3 = res.bootstrap_distribution
print('p-value: ', (  (dlt3 - m >= m) | (dlt3 - m <= -m) ).sum()/M )
#>>p-value:  0.3611361136113611

Этот подход идентичен тому, на который я ссылался в самом начале статьи.

Да, p-значения получились такими, какими мы и ожидали. Но в обоих случаях мы смещаем наше распределение на m. Это потому что мы генерируем наши выборки из двух имитаций генеральных совокупностей. И это, вообще говоря, не такой уж и интерпретируемый результат.

Если верна нулевая гипотеза, то на самом деле у нас одна генеральная совокупность с k1+k2 положительными исходами при n1+n2 испытаниях. Вообще говоря все тесты выше из этого и исходят. Тогда можно модифицировать пример следующим образом:

#Тут нам не нужно ничего никуда смещать.
sales1 = rng.binomial(n1,(k1+k2)/(n1+n2),M)
sales2 = rng.binomial(n2,(k1+k2)/(n1+n2),M)
dlt4 = sales1/n1 - sales2/n2
print('p-value: ', ( (dlt4 >= m) | (dlt4 <= - m) ).sum()/M )
#>>p-value:  0.3625362536253625

px.histogram({'bootstrap':res.bootstrap_distribution-m
              ,'craftbootstrap':dlt4}
             ,title='Распределение разности конверсий'
             ,histnorm='probability'
             ,barmode="overlay"
             ,nbins=500)

Так выглядит наше бутстрап-распределение. Понятно, почему оно изрезанное — т.к. конверсия на самом деле дискретна и это образует вот такую картину.

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

sales1 = rng.binomial(n1,k1/n1,M)
sales2 = rng.binomial(n2,k2/n2,M)
dlt5 = ( sales1/n1 / ( sales2/n2 ) )
px.histogram({'craftbootstrap':dlt5}
             ,title='Распределение отношения конверсий'
             ,histnorm='probability'
             ,barmode="overlay"
             ,nbins=100)

Так же можно получить распределение любой другой функции, если вдруг потребуется

Итого

Бывает смотришь на забор. В нем конечно есть какие-то щелочки, но ничего не видно. Стоит поехать вдоль забора, как сразу вырисовывается картина. Это называется синтез апертуры. Буквально посмотреть на один и тот же объект под разными углами. Цель этой работы как раз в этом и заключается — посмотреть с разных точек зрения на одну и ту же задачу.