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

Вот тут ув. @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 измерения с частотами. Какова вероятность, что распределения частот между ними на самом деле одинаковы, а наблюдаемая разница случайна (принимаем - нулевую гипотезу)? Порешаем прошлый пример.

Rоличество степеней свободы:
В данном случае равно (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 основана на том, что есть несколько категориальных переменных, которые в среднем принимают различные значения. В нашем случае - это тестовая и контрольная группа посетителей на сайт. Есть изменчивость данных, которая складывается из изменчивости внутри групп так и изменчивости между группами (внутригрупповая и межгрупповая дисперсии). На цифрах это вот так выглядит.

Пусть есть группы , где 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 раскладывается на

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

$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$

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

Теперь можно посчитать 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}$ ) и, соответственно, общее количество экспериментов (показов).

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

Первый шаг

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

$\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}$

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

Первый шаг

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

Событие С: К нам пришли 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: . Воспользуемся формулой Байеса и независимостью событий 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)}$

Тут применено соображение что т.к. 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} }$

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

Мы меняем параметр на т.е. в первой группе было не 34 продажи, а 35 или 36 например. При этом должно быть соблюдено условие, где . Т.е. если в первой группе было 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 таких выборок с возвращением в тестовой и контрольной группе. Каждый раз считаем разницу конверсий в группах, т.е. , где . Это такая "логическая имитация" биномиального процесса

#Создали дефолтный генератор
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<k1).sum()
    smpl2 = rng.integers(n2,size=n2)
    k12 = (smpl2<k2).sum()
    dlt[i] = k11/n1 - k12/n2
print('p-value: ', (  (dlt - m >= 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)

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

Итого

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

Источник: https://habr.com/ru/articles/815291/

Вернуться к списку

Интересные статьи

Решите задачу и помогите Саше перестать прокрастинировать

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

Как найти задачу, которую можно решить с помощью машинного обучения

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

Умещаем простую 16-битную VM в 125 строк Си

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

Решаем логистическую задачу: алгоритм привязки фактической и плановой стоянок автомобилей

Привет жителям Хабра! Меня зовут Андрей, я разработчик облачного сервиса по планированию маршрутов Муравьиная логистика. В этой статье я хочу поделиться опытом использования алгоритма, который может п...

Мы разрешаем людям учиться в рабочее время: как мы запустили обучение собственными силами, и что из этого получилось

Привет, Хабр! Меня зовут Лена Хасанова, я — заместитель руководителя в агентстве аутсорс тестирования «Кавычки». Я занимаюсь постановкой процессов, руковожу проектами, ку...