- BrainTools - https://www.braintools.ru -

Байесовские А-Б-тесты: связь с p-значениями

Показана численная близость p-значений t-теста, chi^2-теста и U-критерия Манна-Уитни в А/Б-тестах вероятностям лучшей группы байесовских моделей. Соотношения выполняются несмотря на различия в определениях.

– P-значения [1]
– Т-тест [2]
– Тест chi^2 [3]
– U-критерий Манна-Уитни [4]
– Ссылки [5]

Блокнот: https://github.com/andrewbrdk/Bayesian-AB-Testing/blob/main/appendices/Связь_с_p-значениями.ipynb [6]

import numpy as np
import pandas as pd
import scipy.stats as stats
import plotly.graph_objects as go

np.random.seed(7)

P-значения

P-значения используют в проверках нулевых гипотез. В этом методе формулируют гипотезу H_0 об изучаемом процессе и выбирают статистический тест T – случайную величину с известным распределением P_{T}(x | H_0) в предположении H_0. Считают реализацию величины T в данных – тестовую статистику x_{0}[TestStat [7]]. Вероятность получить фактическое или более экстремальное значение тестовой статистики называют односторонним p-значением p=P_{T}(x ge x_{0} | H_0)[PVal [8], TailedTests [9]]. Если вероятность достаточно мала, гипотезу H_0 отвергают, если нет – оставляют.

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

В предположении гипотезы H_0 распределение тестовой статистики P_{T}(x | H_0). Вероятность получить фактическое x_0 или более экстремальное значение тестовой статистики называют p-значением p=P_{T}(x ge x_{0} | H_0).

Решение о гипотезе H_0 принимают по p-значению p=P_{T}(x ge x_{0} | H_0), тогда как вероятность гипотезы с учетом собранных данных оценивается величиной P(H_0 | x_0). По соотношению Байеса P(H_0 | x_0) propto P_{T}(x=x_{0} | H_0) P(H_0). Т.е. для выбора гипотезы нужно считать вероятности получить данные в рамках конкурирующих гипотез и сравнивать друг с другом с учетом априорных вероятностей.

begin{gather}p=P_{T}(x ge x_{0} | H_0) \ , \ P(H_0 | x_0)=frac{P_{T}(x=x_{0} | H_0) P(H_0)}{P_{T}(x=x_{0} | H_0) P(H_0) + P_{T}(x=x_{0} | {sim}H_0) P({sim}H_0)}end{gather}

В А/Б-тестах распространены t-тест средних, chi^2-тест пропорций и U-критерий Манна-Уитни. Далее показано как при определенных условиях p-значения этих тестов численно близки байесовским вероятностям параметров одной группы больше другой. Соотношения выполняются несмотря на различия в определениях.

Т-тест

Средние сравнивают t-тестами [TTest [10]]. Пусть есть выборки размера N_A, N_B из двух случайных величин A, B. В предположении одинаковых точных средних H_0: E[A]=E[B] для отношения разности выборочных средних к стандартной ошибке [11] разности средних X=overline{Delta}/s_{Delta} ожидают t-распределение [WelchT [12]]. При достаточно большом количестве данных оно близко стандартному нормальному text{Norm}(0, 1) [TDist [13]]. По выборкам считают фактическое отношение x_0=overline{Delta}/s_{Delta}. Вычисляют вероятность получить x_0 или большее значение – одностороннее p-значение P_{X}(x ge x_0 | H_0). Если оно меньше заданного уровня значимости, средние в группах считают неравными.

begin{gather} overline{A}=frac{1}{N_{A}} sum_{i=1}^{N_{A}} A_i, quad s_A^2=frac{1}{N_A} sum_{i=1}^{N_A} (A_i - overline{A})^2, quad text{так же для } B \ X=frac{overline{Delta}}{s_{Delta}}, quad overline{Delta}=overline{B} - overline{A}, quad s^2_{Delta}=frac{s_A^2}{N_A} + frac{s_B^2}{N_B} \ H_0: E[A]=E[B], quad P_{X}(x | H_0) approx text{Norm}(x; 0, 1) \ x_0 - text{реализация } X, quad p=P_{X}(x ge x_0 | H_0) end{gather}

В А/Б-тесте нужно выбрать группу с большим средним. Вместо p-значения P_{X}(x ge x_0 | H_0) интересна вероятность среднего B больше A при условии собранных данных P(mu_B > mu_A | A_i, B_j )=P(mu_{Delta} > 0 | A_i, B_j ), где mu_A, mu_B, mu_{Delta} – оценки точных средних соответствующих распределений. Эту вероятность можно оценить байесовским моделированием. Будет строиться оценка среднего разности mu_{Delta}. По центральной предельной теореме распределение выборочных средних близко нормальному. Поэтому в качестве правдоподобия выбирается нормальное распределение P(overline{Delta} | mu_{Delta})=text{Norm}(overline{Delta} | mu_{Delta}, s_{Delta}^2). Используется модель с одним случайным параметром – средним mu_{Delta}, дисперсия выбирается s_{Delta}^2 на основе данных. Моделирование применяется к выборочным средним, а не исходным выборкам, поэтому для обновления параметров используется только одна точка overline{Delta}. Для сопряженности априорное распределение также выбирается нормальным P(mu_{Delta})=text{Norm}(mu_{Delta} | mu_0, sigma_0^2). Апостериорное распределение будет нормальным с обновленными средним и дисперсией P(mu_{Delta} | overline{Delta})=text{Norm}(mu_{Delta} | mu_{N}, sigma_{N}^2) [ConjPrior [14]]. При достаточно широком априорном распределении sigma_{0}^2 gg s^2_{Delta} апостериорное распределение приближенно P(mu_{Delta} | overline{Delta}) approx text{Norm}(mu_{Delta} | overline{Delta}, s^2_{Delta}). Вероятность среднего в одной группе больше другой можно записать как положительную область нормального распределения с центром в точке x_0=overline{Delta}/s_{Delta} и единичной дисперсией P(mu_B > mu_A | A_i, B_j )=P(mu_{Delta} > 0 | overline{Delta}) approx P(text{Norm}(x > 0 | x_0, 1)).

begin{split} P(overline{Delta} | mu_{Delta}) &=text{Norm}(overline{Delta} | mu_{Delta}, s_{Delta}^2), quad P(mu_{Delta})=text{Norm}(mu_{Delta} | mu_0, sigma_0^2) \ P(mu_{Delta} | overline{Delta}) &=text{Norm}(mu_{Delta} | mu_{N}, sigma_{N}^2), quad sigma_{N}^2=frac{sigma_{0}^2 s_{Delta}^2}{s_{Delta}^2 + sigma_{0}^2}, quad mu_{N}=mu_{0} frac{sigma_{N}^2}{sigma_{0}^2} + frac{sigma_{N}^2}{s_{Delta}^2} overline{Delta} \ mu_0=0, & , sigma_{0}^2 gg s^2_{Delta}: , sigma_N^2 approx s^2_{Delta}, , mu_N approx overline{Delta}, \ P(mu_{Delta} | overline{Delta}) & approx text{Norm}(mu_{Delta} | overline{Delta}, s^2_{Delta}) \ P(mu_B > mu_A | A_i, B_j ) &=P(mu_{Delta} > 0 | overline{Delta}) \ & approx P(text{Norm}(mu_{Delta} > 0 | overline{Delta}, s^2_{Delta}))=P(text{Norm}(x > 0 | x_0, 1)) end{split}

По симметрии [15] нормального распределения P(text{Norm}(x > x_0 | 0, 1))=P(text{Norm}(x < 0 | x_0, 1)). Поэтому p-значение одностороннего t-теста близко вероятности среднего одной группы больше другой.

begin{split} p=P_{X}(x > x_0 | H_0) &=P left( text{Norm}(x > x_0| 0, 1) right) \ &=P left( text{Norm}(x < 0 | x_0, 1) right) \ &=1 - P left( text{Norm}(x > 0 | x_0, 1) right) approx 1 - P(mu_B > mu_A | A_i, B_j ) end{split}

Для демонстрации ниже заданы два нормальных распределения с разными средними. По выборке байесовская оценка вероятности P(mu_B > mu_A | A_i, B_j) сравнивается с p-значением t-теста. Используется односторонний t-тест с разными дисперсиями групп (equal_var=False, alternative) [ScipyTTestInd [16]]. Вероятность
p=P(x > x_0 | H_0)=P left( text{Norm}(x > x_0| 0, 1) right) закрашена темным, P(text{Norm}(x < 0 | x_0, 1)) approx 1 - P(mu_B > mu_A | A_i, B_j) закрашена светлым. По свойствам нормального распределения площади этих областей совпадают. Таким образом p-значение численно близко байесовской оценке вероятности. Стоит помнить, что они не эквивалентны – у них разные определения.

Т-распределение и апостериорное среднего разности 
def posterior_diff_scaled(sampA, sampB, mu0=None, s20=None):
    delta = sampB.mean() - sampA.mean()
    s2delta = sampA.var() / sampA.size + sampB.var() / sampB.size
    mu0 = mu0 or 0
    s20 = s20 or 30 * s2delta
    s2n = s2delta * s20 / (s2delta + s20)
    mun = mu0 * s2n / s20 + delta * s2n / s2delta
    return stats.norm(loc=mun/np.sqrt(s2n), scale=1)

muA = 0.1
muB = 0.115
sigma = 1.3

exactA = stats.norm(muA, sigma)
exactB = stats.norm(muB, sigma)

N = 30000
sampA = exactA.rvs(size=N)
sampB = exactB.rvs(size=N)

a = 'greater' if np.mean(sampA) > np.mean(sampB) else 'less'
t_stat, p_value = stats.ttest_ind(sampA, sampB, equal_var=False, alternative=a)
t_stat = np.abs(t_stat)

post_dist = posterior_diff_scaled(sampA, sampB)
mean_b_gt_a = 1 - post_dist.cdf(0)

xaxis_min = -7
xaxis_max = 8
x = np.linspace(xaxis_min, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=stats.norm.pdf(x, loc=0, scale=1), 
                         line_color='black', opacity=0.8, name='$mathrm{Norm}(0, 1)$'))
fig.add_trace(go.Scatter(x=x[x>t_stat], y=stats.norm.pdf(x[x>t_stat], loc=0, scale=1), 
                         line_color="rgba(0, 0, 0, 0.7)", name='$P(x>x_0 | H_0)$', fill="tozeroy", fillcolor="rgba(0, 0, 0, 0.7)"))
fig.add_trace(go.Scatter(x=x, y=post_dist.pdf(x), 
                         line_color='black', opacity=0.2, name='$mathrm{Norm}(x_0, 1)$'))
fig.add_trace(go.Scatter(x=x[x<0], y=post_dist.pdf(x[x<0]), 
                         line_color="rgba(128, 128, 128, 0.2)", name='$P(mu_{Delta} < 0 | overline{Delta})$', fill="tozeroy", fillcolor="rgba(128, 128, 128, 0.2)"))
fig.add_trace(go.Scatter(x=[0, 0], y=[0, max(stats.norm.pdf(x, loc=0, scale=1))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '0'], textposition="top center", 
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=[t_stat, t_stat], y=[0, max(stats.norm.pdf(x, loc=0, scale=1))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$x_0$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.update_layout(title=r'$Ttext{-распределение и апостериорное среднего разности}$',
                  xaxis_title='$x$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[xaxis_min, xaxis_max],
                  hovermode="x",
                  template="plotly_white",
                  height=500)
fig.show()

print(f'p-value P(x>x0 | H0): {p_value}')
print(f'1 - p: {1 - p_value}')
print(f'Bayes P(pb > pa): {mean_b_gt_a}')
Темная линия - t-распределение, темная закрашенная область - p-значение. Светлая линия - апостериорное распределение среднего разности, светлая закрашенная область - вероятность среднего A меньше B. Площади закрашенных областей совпадают: p-value P(x>x0 | H0): 0.014, 1 - p: 0.986, Bayes P(mu_delta > 0): 0.985

Темная линия – t-распределение, темная закрашенная область – p-значение. Светлая линия – апостериорное распределение среднего разности, светлая закрашенная область – вероятность среднего A меньше B. Площади закрашенных областей совпадают: p-value P(x>x0 | H0): 0.014, 1 – p: 0.986, Bayes P(mu_delta > 0): 0.985

Численную близость p-значения и байесовской вероятности P(mu_{Delta} > 0 | overline{Delta}) можно проверить по количеству правильно угаданных вариантов в серии экспериментов. В каждом эксперименте задается два нормальных распределения. В группе A среднее фиксировано mu=0.1, в B выбирается случайно в диапазоне pm 5% от mu. В группы добавляются данные по n_samp_step точек за шаг. На каждом шаге считается t-тест. Эксперимент останавливается, если 1 - p достигает prob_stop=0.95 или использовано максимальное количество точек n_samp_max. Длительность эксперимента не фиксируется заранее. Миниальный размер выборки n_samp_min + n_samp_step. При остановке эксперимента для сравнения с p-значением считается байесовское апострериорное распределение и вероятность P(mu_{Delta} > 0 | overline{Delta}). Процедура повторяется nexps раз, считается доля правильно угаданных групп во всех экспериментах. Из nexps = 1000 завершено 880, в них правильно угадано 818 вариантов. Точность 0.93 близка prob_stop = 0.95. P-значения в каждом эксперименте близки байесовским вероятностям. При недостаточном n_samp_min доля угаданных вариантов будет меньше prob_stop. Для байесовской модели вместо минимального количества точек можно задать меньшее значение sigma_0^2 – эффект будет близок заданию минимального размера выборки.

Nexp: 1000, Finished: 880, Correct Guesses: 818, Accuracy: 0.930

Доля групп с большим средним, верно угаданных t-тестом 
cmp = pd.DataFrame(columns=['A', 'B', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best_bayes', '1-pval'])

mu = 0.1
nexps = 1000
cmp['A'] = [mu] * nexps
cmp['B'] = mu * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B'] > r['A'] else 'A', axis=1)

n_samp_max = 3_000_000
n_samp_step = 10_000
n_samp_min = 100_000
prob_stop = 0.95

for i in range(nexps):
    muA = cmp.at[i, 'A']
    muB = cmp.at[i, 'B']
    exact_dist_A = stats.norm(loc=muA, scale=1)
    exact_dist_B = stats.norm(loc=muB, scale=1)
    n_samp_current = n_samp_min
    sampA = exact_dist_A.rvs(n_samp_max)
    sampB = exact_dist_B.rvs(n_samp_max)
    post_dist = None
    mean_b_gt_a_bayes = np.nan
    while n_samp_current < n_samp_max:
        n_samp_current += n_samp_step
        a = 'greater' if np.mean(sampA[:n_samp_current]) > np.mean(sampB[:n_samp_current]) else 'less'
        t_stat, p_value = stats.ttest_ind(sampA[:n_samp_current], sampB[:n_samp_current], equal_var=False, alternative=a)
        p_best_t = 1 - p_value
        best_gr = 'A' if p_best_t >= prob_stop and a == 'greater' else 'B' if p_best_t >= prob_stop and a == 'less' else None
        if best_gr:
            post_dist = posterior_diff_scaled(sampA[:n_samp_current], sampB[:n_samp_current])
            mean_b_gt_a_bayes = 1 - post_dist.cdf(0)
            cmp.at[i, 'A_exp'] = sampA[:n_samp_current].mean()
            cmp.at[i, 'B_exp'] = sampB[:n_samp_current].mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_current
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best_bayes'] = max(mean_b_gt_a_bayes, 1 - mean_b_gt_a_bayes)
            cmp.at[i, '1-pval'] = 1 - p_value
            break
    print(f'done {i}: nsamp {n_samp_current}, best_gr {best_gr}, Bayes P(b>a) {mean_b_gt_a_bayes:.4f}, T-test p-val {p_value:.4f}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(30))
finished = np.sum(cmp['best_exp'].notna())
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Finished: {finished}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / finished}")
Байесовские А-Б-тесты: связь с p-значениями - 85

Тест Хи-квадрат

Конверсии могут сравнивать chi^2-тестом [Chi2Test [17]]. Статистика chi^2 Пирсона [Chi2Pearson [18]] для мультиномиальных распределений определена chi^2=sum_{i=1}^k (S_i - Np_i)^2/Np_i, где N – общее количество наблюдений, S_i и N p_i – фактическое и ожидаемое количество наблюдений i-категории при доле i-категории p_i. Для биномиального распределения chi^2=(S - Np)^2/N p (1-p). По центральной предельной теореме (S - Np)/sqrt{N p (1-p)} стремится к стандартному нормальному распределению, квадрат этой величины совпадает с chi^2. Распределение суммы квадратов k нормальных случайных величин называют chi^2-распределением с k степенями свободы chi^2_k=sum_{i=1}^{k} X_i^2,, X_i sim text{Norm}(0,1) [Chi2Dist [19]]. Поэтому статистика chi^2 стремится к chi_1^2-распределению.

begin{split} chi^2 &=sum_{i=1}^k frac{(S_i - Np_i)^2}{N p_i} \ &=frac{(S - N p)^2}{N p} + frac{((N - S) - N (1-p))^2}{N (1-p)} \ &=frac{(S - Np)^2}{N p (1-p)} to chi_1^2, quad N to infty \ chi^2_k &=sum_{i=1}^{k} X_i^2,, X_i sim text{Norm}(0,1) end{split}

Для А/Б-теста конверсий с двумя группами в предположении одинаковых ожидаемых конверсий p=(S_A + S_B)/(N_A + N_B) тестовую статистику можно привести к виду chi^2=Delta p^2/s_{Delta}^2, Delta p=S_B/N_B - S_A/N_A, s_{Delta}^2=p(1-p) (1 / N_A + 1 / N_B). При большом количестве точек распределение chi^2 близко chi_1^2. Т.к. распределение chi^2_1 получается возведением в квадрат стандартного нормального распределения, область p-значения p=P_{chi_1^2}(x > chi^2 | H_0) соответствует областям Pleft(text{Norm}(x > chi cup x < -chi | 0, 1)right). По симметрии нормального распределения площади x > chi и x < -chi одинаковы.

begin{gather} A sim mathrm{Bernoulli}(p_A), , S_A=sum_{i=1}^{N_A} A_i, quad B sim mathrm{Bernoulli}(p_B), , S_B=sum_{i=1}^{N_B} B_i, \ p=frac{S_A + S_B}{N_A + N_B}, quad Delta p=frac{S_B}{N_B} - frac{S_A}{N_A}, quad s_{Delta}^2=frac{p(1-p)}{N_A} + frac{p(1-p)}{N_B} \ begin{split} H_0: E[A]=E[B], quad chi^2 &=frac{(S_A - N_A p)^2}{N_A p (1-p)} + frac{(S_B - N_B p)^2}{N_B p (1-p)} \ &=frac{N_A N_B Delta p^2}{(N_A + N_B) p (1-p)} \ &=frac{Delta p^2}{s_{Delta}^2} to chi_1^2, , n to infty end{split} \ , \ begin{split} text{p-val} &=P_{chi_1^2}(x > chi^2 | H_0) \ &=P left( text{Norm}(x > chi cup x < -chi; 0, 1) right) \ &=2 Pleft( text{Norm}(x > chi; 0, 1) right) end{split} end{gather}

Для байесовской оценки вероятности конверсии одной группы больше другой P(theta_B > theta_A | A_i, B_j) правдоподобие в каждой группе задается биномиальным распределением P(S | theta, N)=mbox{Binom}(S | theta, N), априорное – бета-распределением P(theta)=mbox{Beta}(theta; alpha, beta). Апостериорное будет бета-распределением с обовленными параметрами P(theta | S, N)=mbox{Beta}(theta; alpha + S, beta + N - S). При типичных значениях N, S бета-распределение близко нормальному. Распределение разности конверсий также будет нормальным. При равномерных априорных распределениях и слабой разнице между группами P_{Delta theta}(x) approx mbox{Norm}left(x; Delta p, s_{Delta}^2 right) с параметрами Delta p, s_{Delta}^2 как в chi^2-тесте. Вероятность конверсии одной группы больше другой будет площадью этого распределения в положительной области P(theta_B > theta_A | A_i, B_j)=P(Delta theta > 0 | A_i, B_j) approx Pleft(text{Norm}(x > 0; chi, 1) right).

begin{split} P(S | theta, N) &=mbox{Binom}(S | theta, N) \ P(theta) &=mbox{Beta}(theta; alpha, beta) \ P(theta | S, N) &=mbox{Beta}(theta; alpha + S, beta + N - S) \ & approx mbox{Norm}(theta; mu, sigma^2), quad mu=S / N, quad sigma^2=mu (1 - mu) / N, quad S, N gg alpha, beta gg 1 end{split}begin{split} P_{theta_A}(x) &=mbox{Beta}(x; alpha_A + S_A, beta_A + N_A - S_A) approx mbox{Norm}(x; mu_A, sigma^2_A), \ P_{theta_B}(x) &=mbox{Beta}(x; alpha_B + S_B, beta_B + N_B - S_B) approx mbox{Norm}(x; mu_B, sigma^2_B) \ P_{Delta theta}(x) & approx mbox{Norm}left(x; mu_B - mu_A, sigma_A^2 + sigma_B^2right) approx mbox{Norm}left(x; Delta p, s_{Delta}^2 right), quad text{при } p approx S_A/N_A approx S_B/N_B end{split}begin{split}P(theta_B > theta_A | A_i, B_j) &=P(Delta theta > 0 | A_i, B_j) \ & approx Pleft(text{Norm}(x > 0; Delta p, s_{Delta}^2) right)=Pleft(text{Norm}(x > 0; chi, 1) right), quad chi=Delta p / s_{Delta} end{split}

По симметрии Pleft( text{Norm}(x > chi | 0, 1) right)=Pleft(text{Norm}(x < 0; chi, 1) right). Поэтому P(theta_B > theta_A | A_i, B_j) approx 1 - mbox{p-val}/2.

begin{gather} begin{split} text{p-val} &=P_{chi_1^2}(x > chi^2 | H_0) \ &=2 Pleft( text{Norm}(x > chi; 0, 1) right) \ &=2 Pleft(text{Norm}(x < 0; chi, 1) right) \ & approx 2 left( 1 - P(theta_B > theta_A | A_i, B_j) right) end{split} \ , \ P(theta_B > theta_A | A_i, B_j) approx 1 - frac{text{p-val}}{2} end{gather}

Соотношение P(theta_B > theta_A | A_i, B_j) approx 1 - mbox{p-val}/2 проверяется по выборке из двух распределений Бернулли с конверсиями p_A=0.1 и p_B=0.103. Данные для chi^2-теста задаются в виде таблицы со строками S_A, N_A-S_A и S_B, N_B - S_B [ScipyChi2Con [20]]. P-значение chi^2-теста не позволяет выбрать между p_A > p_B и p_B > p_A, поэтому дополнительно сравниваются средние в выборках. Распределение chi^2_1 на первом графике ниже. Закрашенная область соответствует p-значению mbox{p-val}=P_{chi_1^2}(x > chi^2 | H_0). На втором графике закрашенные темные области x > chi и x < - chi соответствуют области p-значения при возведении в квадрат. Серый график – нормальное распределение text{Norm}(chi, 1). Закрашенная область серого графика приближенно равна 1 - P(theta_B > theta_A | A_i, B_j). Площади закрашенной серой и каждой из темных областей совпадают. Связь p-значения и байесовской оценки вероятности численно выполняется.

Распределение Хи-квадрат, p-значение и байесовская модель 
#np.random.seed(27)

def posterior_dist_binom(ns, ntotal, a_prior=1, b_prior=1):
    a = a_prior + ns
    b = b_prior + ntotal - ns 
    return stats.beta(a=a, b=b)

def diff_ba_scaled(post_dist_A, post_dist_B):
    m = post_dist_B.mean() - post_dist_A.mean()
    v = post_dist_A.var() + post_dist_B.var()
    return stats.norm(loc=m/np.sqrt(v), scale=1)


pA = 0.1
pB = pA * 1.03

exactA = stats.bernoulli(pA)
exactB = stats.bernoulli(pB)

N = 10000
sampA = exactA.rvs(size=N)
sampB = exactB.rvs(size=N)
SA = np.sum(sampA)
SB = np.sum(sampB)

t = np.array([
    [SA,     N - SA],
    [SB,     N - SB]
])
chi2_stat, p_value_chi2, dof, expected = stats.chi2_contingency(t, correction=False)
chi = np.sqrt(chi2_stat)
p_A_samp = SA / N
p_B_samp = SB / N
pb_gt_pa_chi = 1 - p_value_chi2 / 2
pb_gt_pa_chi = pb_gt_pa_chi if p_B_samp > p_A_samp  else 1 - pb_gt_pa_chi

post_dist_A = posterior_dist_binom(ns=SA, ntotal=N)
post_dist_B = posterior_dist_binom(ns=SB, ntotal=N)
diff_dist_scaled = diff_ba_scaled(post_dist_A, post_dist_B)
pb_gt_pa_bayes = 1 - diff_dist_scaled.cdf(0)

xaxis_min = 0
xaxis_max = 5
x = np.linspace(xaxis_min, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=stats.chi.pdf(x, df=1), 
                         line_color='black', opacity=0.8, name=f'$chi^2_1$'))
fig.add_trace(go.Scatter(x=[chi2_stat, chi2_stat], y=[0, max(stats.chi.pdf(x, df=1))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$chi^2$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=x[x>chi2_stat], y=stats.chi.pdf(x[x>chi2_stat], df=1), 
                         line_color='black', opacity=0.8, name='$P_{chi_1^2}(x > chi^2)$', fill="tozeroy", fillcolor="rgba(0, 0, 0, 0.7)"))
fig.update_layout(title='Хи-квадрат',
                  xaxis_title='$x$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[xaxis_min, xaxis_max],
                  hovermode="x",
                  template="plotly_white",
                  height=500)
fig.show()

xaxis_min = -5
xaxis_max = 5
x = np.linspace(xaxis_min, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=stats.norm.pdf(x, loc=0, scale=1), 
                         line_color='black', opacity=0.8, name='$mathrm{Norm}(0, 1)$'))
fig.add_trace(go.Scatter(x=[0, 0], y=[0, max(stats.norm.pdf(x, loc=0, scale=1))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '0'], textposition="top center", 
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=[chi, chi], y=[0, max(stats.norm.pdf(x, loc=0, scale=1))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$chi$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=[-chi, -chi], y=[0, max(stats.norm.pdf(x, loc=0, scale=1))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$-chi$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=x[x>chi], y=stats.norm.pdf(x[x>chi], loc=0, scale=1), 
                         line_color='black', opacity=0.8, name='$P(mathrm{Norm}(x > chi cup x < -chi | 0, 1)$', fill="tozeroy", fillcolor="rgba(0, 0, 0, 0.7)"))
fig.add_trace(go.Scatter(x=x[x<-chi], y=stats.norm.pdf(x[x<-chi], loc=0, scale=1), 
                         line_color='black', opacity=0.8, name='$P(mathrm{Norm}(x > chi | 0, 1)$', fill="tozeroy", fillcolor="rgba(0, 0, 0, 0.7)",
                         showlegend=False))
fig.add_trace(go.Scatter(x=x, y=diff_dist_scaled.pdf(x), 
                         line_color='black', opacity=0.2, name=r'$mathrm{Norm}(chi, 1)$'))
fig.add_trace(go.Scatter(x=x[x<0], y=diff_dist_scaled.pdf(x[x<0]), 
                         line_color="rgba(128, 128, 128, 0.2)", name=r'$P(mathrm{Norm}(x < 0 | chi, 1))$', fill="tozeroy", fillcolor="rgba(128, 128, 128, 0.2)"))
fig.update_layout(title=r'$Ptext{-значение и байесовская модель}$',
                  xaxis_title='$x$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[xaxis_min, xaxis_max],
                  hovermode="x",
                  template="plotly_white",
                  height=500)
fig.show()

print(f'Bayes P(pb > pa): {pb_gt_pa_bayes:.5g}')
print(f"Chi2: 1-pval/2:   {pb_gt_pa_chi:.5g}")
Распределение , закрашенная область соответствует p-значению.

Распределение chi^2, закрашенная область соответствует p-значению.
Темная линия - стандартное нормальное распределение, при возведении в квадрат получится -распределение. Темные закрашенные области соответствуют p-значению -распределения. Светлая линия - апостериорное байесовское распределение разности конверсий. Серая область - вероятность конверсии группы А больше Б. Площадь серой закрашенной области совпадает с темными закрашенными областями. Bayes P(theta_b > theta_a): 0.916, Chi2: 1-pval/2: 0.916.

Темная линия – стандартное нормальное распределение, при возведении в квадрат получится chi^2-распределение. Темные закрашенные области соответствуют p-значению chi^2-распределения. Светлая линия – апостериорное байесовское распределение разности конверсий. Серая область – вероятность конверсии группы А больше Б. Площадь серой закрашенной области совпадает с темными закрашенными областями. Bayes P(theta_b > theta_a): 0.916, Chi2: 1-pval/2: 0.916.

P-значение chi^2-теста используется для выбора вариантов с большей конверсией в серии экспериментов. В каждом эксперименте 2 группы, конверсия p_A=0.1 фиксирована, p_B выбирается случайно в диапазоне pm5% от p_A. В каждой группе добавляются данные по n_samp_step точек за шаг. На каждом шаге вычисляется p-значение chi^2-теста и P(theta_B > theta_A | A_i, B_j). Эксперимент останавливается, если оценка вероятности конверсии одной группы больше другой превышает prob_stop или набрано максимальное количество точек n_samp_max. Всего из nexps=1000 завершено 986, верно угадано 940 варианта. Доля 0.953 близка prob_stop = 0.95.

Nexp: 1000, Finished: 986, Correct Guesses: 940, Accuracy: 0.953

Доля групп с большей конверсией, верно угаданных chi2-тестом 
cmp = pd.DataFrame(columns=['A', 'B', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best_bayes', 'p_best_chi'])

p = 0.1
nexps = 1000
cmp['A'] = [p] * nexps
cmp['B'] = p * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B'] > r['A'] else 'A', axis=1)

n_samp_max = 5_000_000
n_samp_step = 10_000
prob_stop = 0.95

for i in range(nexps):
    pA = cmp.at[i, 'A']
    pB = cmp.at[i, 'B']
    exact_dist_A = stats.bernoulli(p=pA)
    exact_dist_B = stats.bernoulli(p=pB)
    n_samp_total = 0
    ns_A = 0
    ns_B = 0
    while n_samp_total < n_samp_max:
        dA = exact_dist_A.rvs(n_samp_step)
        dB = exact_dist_B.rvs(n_samp_step)
        n_samp_total += n_samp_step
        ns_A = ns_A + np.sum(dA)
        ns_B = ns_B + np.sum(dB)
        p_A_samp = ns_A / n_samp_total
        p_B_samp = ns_B / n_samp_total
        t = np.array([
            [ns_A,     n_samp_total - ns_A],
            [ns_B,     n_samp_total - ns_B]
        ])
        chi2_stat, p_value_chi, dof, expected = stats.chi2_contingency(t, correction=False)
        pb_gt_pa_chi = 1 - p_value_chi / 2
        pb_gt_pa_chi = pb_gt_pa_chi if p_B_samp > p_A_samp  else 1 - pb_gt_pa_chi
        best_gr = 'B' if pb_gt_pa_chi >= prob_stop else 'A' if 1 - pb_gt_pa_chi >= prob_stop else None
        if best_gr:
            post_dist_A = posterior_dist_binom(ns=ns_A, ntotal=n_samp_total)
            post_dist_B = posterior_dist_binom(ns=ns_B, ntotal=n_samp_total)
            pb_gt_pa_bayes = 1 - diff_ba_scaled(post_dist_A, post_dist_B).cdf(0)
            cmp.at[i, 'A_exp'] = p_A_samp
            cmp.at[i, 'B_exp'] = p_B_samp
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best_bayes'] = max(pb_gt_pa_bayes, 1 - pb_gt_pa_bayes)
            cmp.at[i, 'p_best_chi'] = max(pb_gt_pa_chi, 1 - pb_gt_pa_chi)
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P_best Bayes {max(pb_gt_pa_bayes, 1 - pb_gt_pa_bayes):.4f}, Chi (1-pval/2): {1 - p_value_chi/2:.4f}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(30))
finished = np.sum(cmp['best_exp'].notna())
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Finished: {finished}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / finished}")
Байесовские А-Б-тесты: связь с p-значениями - 167

U-критерий Манна-Уитни

Для выборок размера N_A, N_B из двух случайных величин A, B статистика Манна-Уитни [MannWhitneyU [21]] определена попарным сравнением элементов. Для непрерывных распределений вероятность совпадения элементов в выборках нулевая. В этом случае U-статистика определена как количество пар (A_i, B_j), где элемент A_i больше B_j: U_A=sum_{i=1}^{N_A} sum_{j=1}^{N_B} I(A_i > B_j), I – индикаторная функция. Статистику также можно записать в виде U_A=R_A - N_A (N_A + 1)/2, где R_A – сумма рангов элементов A в объединенной выборке. Эквивалентность определений можно увидеть следующим образом: слагаемое N_A (N_A + 1)/2 соответствует минимальной сумме рангов если все элементы A_i меньше B_j и считается как сумма арифметической прогрессии. Если наибольший элемент A_i больше n элементов B_j, то U_A=n и R_A=N_A (N_A + 1)/2 + n. При большом количестве точек отношение U_A к общему количеству пар U_A / N_A N_B стремится к вероятности точки из распределения A больше B.

begin{split} A, B & - text{непрерывные распределения} \ U_A &=sum_{i=1}^{N_A} sum_{j=1}^{N_B} I(A_i > B_j), quad I(cdot)=1 text{ если условие выполнено, иначе } 0 \ U_A &=R_A - N_A (N_A + 1)/2, quad R_A text{- сумма рангов элементов A в объединенной выборке } A_i, B_j \ frac{U_A}{N_A N_B} & to P(A > B), quad N_A, N_B to infty end{split}

В предположении одинаковых распределений A=B можно посчитать среднее E[U_A] и дисперсию text{Var}(U_A) U-статистики. Величина (U_A - E[U_A])/sqrt{text{Var}(U_A)} будет стремиться к нормальному распределению. От количества пар удобно перейти к вероятности u_A=U_A / N_A N_B. P-значение вычисляется как p=P left( mbox{Norm}left(x < u_A; 1/2, text{Var}(U_A)/(N_A N_B)^2 right) right) при u_A < 0.5 или как площадь области x > u_A при u_A > 0.5.

begin{gather} H_0: A=B, quad E[U_A]=frac{N_A N_B}{2}, quad text{Var}(U_A)=frac{N_A N_B (N_A + N_B + 1)}{12} \ frac{U_A - E[U_A]}{sqrt{text{Var}(U_A)}} to text{Norm}(0,1), quad N_A, N_B to infty \ u_A=frac{U_A}{N_A N_B}, quad frac{u_A - 1/2}{sqrt{text{Var}(U_A)/(N_A N_B)^2}} to text{Norm}(0,1), quad N_A, N_B to infty \ p=P( x > u_A | H_0 )=P left( mbox{Norm}left(x < u_A; 1/2, text{Var}(U_A)/(N_A N_B)^2 right) right), quad u_A < 0.5 end{gather}

В байесовском подходе вероятность P(B>A) можно оценить сравнением апостериорных предиктивных распределений. Для этого нужно задать модели исходных распределений, построить апостериорные распределения параметров, апостериорные предиктивные распределения и по ним оценить P(B>A). Другой вариант, не требующий предположений о распределениях – моделировать вероятность theta=P(B > A) точки B больше A. Данными будут пары (A_i, B_j), в каждой паре точка B_j больше A_i с вероятностью theta, по общему количеству пар N и количеству пар S с B_j больше A_i нужно оценить theta. Такая задача похожа на моделирование конверсий. Если при формировании пар из исходных выборок каждую точку A сравнивать с каждой точкой B, пары с одинаковыми A_i или B_j будут зависимы. Для независимости каждая точка A_i сравнивается только с одной точкой B_j. Аналогично конверсиям правдоподобие задается биномиальным распределением P(S | theta)=text{Binom}(S | theta, N), априорное – бета-распределением. Апостериорное также будет бета-распределением. Вероятность точки одной группы больше другой равна области theta больше 0.5: P(B > A) approx P(theta > 0.5 | S, N).

begin{gather} theta=P(B > A), quad S=sum_{i=1}^N I(B_i > A_i), quad N_A=N_B=N \ begin{split} P(S | theta) &=text{Binom}(S; theta, N) \ P(theta) &=text{Beta}(theta; alpha, beta) \ P(theta | S) &=text{Beta}(theta; S + alpha, N - S + beta) \ & approx text{Norm}(theta; mu, sigma^2), quad mu=(S + alpha ) / (N + alpha + beta), quad sigma^2=mu (1 - mu) / N, quad N, S gg alpha, beta gg 1 end{split} \ , \ begin{split} P(B > A) approx P(theta > 0.5 | S, N) &=P(text{Norm}(x > 0.5; mu, sigma^2)) end{split} end{gather}

В U-статистике каждая точка A_i сравнивается с каждой B_j, в байесовской модели точки A_i, B_j участвуют только в одном сравнении. Из-за меньшего использования данных дисперсия байесовской модели шире U-статистики. P-значение не будет близко совпадать с байесовской вероятностью P(theta > 0.5 | S, N). При большом количестве точек text{Var}(U_A)/(N_A N_B)^2 to 1/6N, N_A=N_B=N to infty. В байесовской модели при слабом отличии между группами mu approx 0.5, sigma^2 approx 1/4N. Поэтому p lesssim 1 - P(B > A).

begin{split} p & approx P left( mbox{Norm}left(x < u_A; 0.5, 1/6N right) right) \ & lesssim P(text{Norm}(x < 0.5; mu, 1/4N)) \ &=P(A > B) \ & approx 1 - P(B > A) end{split}

Проверка связи p-значения с байесовской вероятностью сделана на примере нормальных распределений. Вероятность точки из распределения B больше A равна вероятности разности B-A больше нуля P(B > A)=P(B - A > 0). Разность нормальных распределений также нормальное распределение со средним, равным разности исходных средних, и дисперсией, равной сумме дисперсий [NormalSum [22]].

begin{gather} A sim text{Norm}(mu_A, sigma_A^2), quad B sim text{Norm}(mu_B, sigma_B^2) \ B - A sim text{Norm}(mu_B - mu_A, sigma_A^2 + sigma_B^2) \ P(B > A)=P(B - A > 0)=1 - F_{B-A}(0) end{gather}

На первом графике ниже показаны исходные нормальные распределения со средними mu1=0, mu2=0.1 и единичной дисперсией. Темная линия на втором графике – U/N_A N_B в предположении эквивалентности распределений. Темная закрашенная область соответствует p-значению [ScipyMannWhitneyU [23]]. Светлая линия показывает байесовскую модель, светлая закрашенная область – байесовскую вероятность P(A>B). Байесовская модель шире U/N_A N_B из-за большей дисперсии. Площадь светлой области обычно больше темной, что соответствует p lesssim P(A>B).

U-статистика и байесовская модель 
#np.random.seed(34)

mu1, sigma1 = 0.0, 1
mu2, sigma2 = 0.1, 1
exactA = stats.norm(loc=mu1, scale=sigma1)
exactB = stats.norm(loc=mu2, scale=sigma2)
p_b_gt_a_exact = 1 - stats.norm.cdf(0, loc=mu2-mu1, scale=np.sqrt(sigma1**2 + sigma2**2))

nA = nB = 1000
sampA = exactA.rvs(nA)
sampB = exactB.rvs(nB)

U, pval = stats.mannwhitneyu(sampA, sampB, alternative='greater')
pval = pval if U / (nA * nB) > 0.5 else 1 - pval
varu = nA * nB * (nA + nB + 1) / 12
p_u = stats.norm(loc=0.5, scale=np.sqrt(varu/(nA*nA*nB*nB)))
ua = U / (nA*nB)

a0 = 1
b0 = 1
S = np.sum(sampB > sampA)
post_theta = stats.beta(a0 + S, b0 + nB - S)


x = np.linspace(-7, 7, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(
    x=x, y=exactA.pdf(x),
    mode='lines', name='A', line_color='black'))
fig.add_trace(go.Scatter(
    x=x, y=exactB.pdf(x),
    mode='lines', name='B', line_color='blue', opacity=0.7))
fig.update_layout(
    title="A, B",
    template="plotly_white"
)
fig.show()

xaxis_min = 0.4
xaxis_max = 0.6
x = np.linspace(xaxis_min, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=p_u.pdf(x), 
                         line_color='black', opacity=0.8, name=f'$U_A/N_A N_B$'))
fig.add_trace(go.Scatter(x=x[x<ua], y=p_u.pdf(x[x<ua]), 
                         line_color='black', opacity=0.8, name=r'$p = P(x < u_A | H_0)$', 
                         fill="tozeroy", fillcolor="rgba(0, 0, 0, 0.7)"))
fig.add_trace(go.Scatter(x=x, y=post_theta.pdf(x),
                         line_color='black', opacity=0.3, name=r'$P(theta | S, N)$'))
fig.add_trace(go.Scatter(x=x[x<0.5], y=post_theta.pdf(x[x<0.5]), 
                         name=r'$P(theta < 0.5 | S, N)$',
                         line_color="rgba(128, 128, 128, 0.3)",
                         fill="tozeroy", fillcolor="rgba(0, 0, 0, 0.3)"))
fig.add_trace(go.Scatter(x=[1-ua, 1-ua], y=[0, max(p_u.pdf(x))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$1-u_A$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=[ua, ua], y=[0, max(p_u.pdf(x))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$u_A$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.add_trace(go.Scatter(x=[0.5, 0.5], y=[0, max(p_u.pdf(x))*1.1], 
                         line_color='black', 
                         mode='lines+text', text=['', '$0.5$'], textposition="top center",
                         line_dash='dash', showlegend=False))
fig.update_layout(
    title=r'$Utext{-статистика и байесовская модель}$',
    template="plotly_white"
)
fig.show()

print(f'CDF(u_a): {p_u.cdf(ua):.5f}')
print(f'scipy min(pval, 1-pval): {min(pval, 1-pval):.5f}')
print(f'Bayes P(theta < 0.5) {post_theta.cdf(0.5):.5f}')
print()
print(f'P(B>A) exact: {p_b_gt_a_exact:.5f}')
print(f'1 - U/(NA NB): {1-ua:.5f}')
print(f'Bayes E[theta] {post_theta.mean():.5f}')
Исходные распределения

Исходные распределения
Темная линия соответствует распределению , закрашенная темная область - p-значению U-теста. Светлая линия - байесовская модель, закрашенная светлая область - оценка вероятности точки из А больше Б. При большом количестве точек p-значение меньше байесовской вероятности описанной модели. pval: 0.011, Bayes P(theta < 0.5) 0.015.

Темная линия соответствует распределению U_A/N_A N_B, закрашенная темная область – p-значению U-теста. Светлая линия – байесовская модель, закрашенная светлая область – оценка вероятности точки из А больше Б. При большом количестве точек p-значение меньше байесовской вероятности описанной модели. pval: 0.011, Bayes P(theta < 0.5) 0.015.

Для байесовской модели и U-статистики проверяется доля правильно угаданных вариантов с большей вероятностью P(B>A) в серии экспериментов. В каждом эксперименте задается два нормальных распределения. Среднее A фиксировано mua = 0.1, среднее B выбирается случайно в пределах pm 5% от A. Вначале проверяется байесовская модель. Всего проводится nexp экспериментов, в группы добавляются данные по n_samp_step точек за раз. На каждом шаге считается байесовсая вероятность P(B > A)=P(theta > 0.5 | S). Эксперимент останавливается, если P(B > A) или P(A > B) превышает prob_stop. В обоих подходах – байесовском и U-тесте – нужно задавать минимальный размер выборки для выполнения заданного уровня точности. В байесовской модели вместо минимальной выборки задаются априорные параметры a0 = 100000, b0 = 100000. Из 1000 экспериментов завершено 860, корректно угадано 827. Точность 0.96 близка заданной prob_stop = 0.95. При остановке в каждом эксперименте считается U-статистика, она как правило превышает байесовскую вероятность.

Nexp: 1000, Finished: 860, Correct Guesses: 827, Accuracy: 0.962

Доля групп, верно угаданных байесовской моделью 
cmp = pd.DataFrame(columns=['A', 'B', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best_bayes', 'p_best_u'])

mua = 0.1
nexps = 1000
cmp['A'] = [mua] * nexps
cmp['B'] = mua * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B'] > r['A'] else 'A', axis=1)

n_samp_max = 5_000_000
n_samp_step = 10000
prob_stop = 0.95

for i in range(nexps):
    mua = cmp.at[i, 'A']
    mub = cmp.at[i, 'B']
    exact_dist_A = stats.norm(loc=mua)
    exact_dist_B = stats.norm(loc=mub)
    n_samp_total = 0
    sampA = exact_dist_A.rvs(n_samp_max)
    sampB = exact_dist_B.rvs(n_samp_max)
    while n_samp_total < n_samp_max:
        n_samp_total += n_samp_step
        a0 = 100000
        b0 = 100000
        Ub = np.sum(sampB[:n_samp_total] > sampA[:n_samp_total])
        post_u_ewise = stats.beta(a0 + Ub, b0 + n_samp_total - Ub)
        pb_gt_pa_bayes = 1 - post_u_ewise.cdf(0.5)
        best_gr = 'B' if pb_gt_pa_bayes >= prob_stop else 'A' if 1 - pb_gt_pa_bayes >= prob_stop else None
        if best_gr:
            U, pval = stats.mannwhitneyu(sampA[:n_samp_total], sampB[:n_samp_total], alternative='greater')
            cmp.at[i, 'A_exp'] = sampA[:n_samp_total].mean()
            cmp.at[i, 'B_exp'] = sampB[:n_samp_total].mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best_bayes'] = max(pb_gt_pa_bayes, 1 - pb_gt_pa_bayes)
            cmp.at[i, 'p_best_u'] = max(pval, 1 - pval)
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P_best Bayes {pb_gt_pa_bayes:.4f}, U p-val: {pval:.4f}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(30))
finished = np.sum(cmp['best_exp'].notna())
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Finished: {finished}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / finished}")
Байесовские А-Б-тесты: связь с p-значениями - 265

Для проверки правильно угаданных вариантов по U-статистике выставлен больший шаг – добавляется по n_samp_step = 200_000 точек в каждой группе. Такой шаг также учитывает минимальный размер выборки. Из 300 экспериментов завершено 262, верно угадано 249 групп. Доля верно угаданных вариантов 0.95 совпала с prob_stop=0.95.

Nexp: 300, Finished: 262, Correct Guesses: 249, Accuracy: 0.950

Доля групп, верно угаданных U-тестом 
cmp = pd.DataFrame(columns=['A', 'B', 'best_exact', 'exp_samp_size', 'A_exp', 'B_exp', 'best_exp', 'p_best_bayes', 'p_best_u'])

mua = 0.1
nexps = 300
cmp['A'] = [mua] * nexps
cmp['B'] = mua * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))
cmp['best_exact'] = cmp.apply(lambda r: 'B' if r['B'] > r['A'] else 'A', axis=1)

n_samp_max = 5_000_000
n_samp_step = 200_000
prob_stop = 0.95

for i in range(nexps):
    mua = cmp.at[i, 'A']
    mub = cmp.at[i, 'B']
    exact_dist_A = stats.norm(loc=mua)
    exact_dist_B = stats.norm(loc=mub)
    n_samp_total = 0
    sampA = exact_dist_A.rvs(n_samp_max)
    sampB = exact_dist_B.rvs(n_samp_max)
    while n_samp_total < n_samp_max:
        n_samp_total += n_samp_step
        U, pval = stats.mannwhitneyu(sampA[:n_samp_total], sampB[:n_samp_total], alternative='greater')
        pb_gt_pa_u = 1 - U / n_samp_total / n_samp_total
        best_gr = 'B' if pval >= prob_stop else 'A' if 1 - pval >= prob_stop else None        
        if best_gr:
            a0 = 100000
            b0 = 100000
            Ub = np.sum(sampB[:n_samp_total] > sampA[:n_samp_total])
            post_u_ewise = stats.beta(a0 + Ub, b0 + n_samp_total - Ub)
            pb_gt_pa_bayes = 1 - post_u_ewise.cdf(0.5)
            cmp.at[i, 'A_exp'] = sampA[:n_samp_total].mean()
            cmp.at[i, 'B_exp'] = sampB[:n_samp_total].mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best_bayes'] = max(pb_gt_pa_bayes, 1 - pb_gt_pa_bayes)
            cmp.at[i, 'p_best_u'] = max(pval, 1 - pval)
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P_best Bayes {pb_gt_pa_bayes:.4f}, U p-val: {pval:.4f}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(30))
finished = np.sum(cmp['best_exp'].notna())
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Finished: {finished}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / finished}")
Байесовские А-Б-тесты: связь с p-значениями - 267

Для t-теста, chi^2-теста и U-критерия Манна-Уитни показаны байесовские модели с вероятностью лучшей группы численно близкой p-значениям. Соотношения выполняются несмотря на различия в определениях p-значений и байесовских вероятностей.

Ссылки

[Chi2Dist] – Chi-squared Distribution [19], Wikipedia.
[Chi2Pearson] – Pearson’s Chi-squared Test [18], Wikipedia.
[Chi2Test] – Chi-squared Test [17], Wikipedia.
[ConjPrior] – Conjugate Prior [14], Wikipedia.
[MannWhitneyU] – Mann–Whitney U Test [21], Wikipedia.
[NormalSum] – Sum of Normally Distributed Random Variables [22], Wikipedia.
[PVal] – P-value [8], Wikipedia.
[ScipyChi2Con] – scipy.stats.chi2_contingency [20], SciPy Reference.
[ScipyMannWhitneyU] – scipy.stats.mannwhitneyu [23], SciPy Reference.
[ScipyTTestInd] – scipy.stats.ttest_ind [16], SciPy Reference.
[TailedTests] – One- and Two-tailed Tests [9], Wikipedia.
[TDist] – Student’s t-distribution [13], Wikipedia.
[TestStat] – Test Statistic [7], Wikipedia.
[TTest] – Student’s t-test [10], Wikipedia.
[WelchT] – Welch’s t-test [12], Wikipedia.

Автор: andrew_brdk

Источник [24]

Сайт-источник BrainTools: https://www.braintools.ru

Путь до страницы источника: https://www.braintools.ru/article/24629

URLs in this post:

[1] P-значения: #pval

[2] Т-тест: #ttest

[3] Image: #chitest

[4] U-критерий Манна-Уитни: #utest

[5] Ссылки: #refs

[6] https://github.com/andrewbrdk/Bayesian-AB-Testing/blob/main/appendices/Связь_с_p-значениями.ipynb: https://github.com/andrewbrdk/Bayesian-AB-Testing/blob/main/appendices/%D0%A1%D0%B2%D1%8F%D0%B7%D1%8C_%D1%81_p-%D0%B7%D0%BD%D0%B0%D1%87%D0%B5%D0%BD%D0%B8%D1%8F%D0%BC%D0%B8.ipynb

[7] TestStat: https://en.wikipedia.org/wiki/Test_statistic

[8] PVal: https://en.wikipedia.org/wiki/P-value

[9] TailedTests: https://en.wikipedia.org/wiki/One-_and_two-tailed_tests

[10] TTest: https://en.wikipedia.org/wiki/Student%27s_t-test

[11] ошибке: http://www.braintools.ru/article/4192

[12] WelchT: https://en.wikipedia.org/wiki/Welch%27s_t-test

[13] TDist: https://en.wikipedia.org/wiki/Student%27s_t-distribution

[14] ConjPrior: https://en.wikipedia.org/wiki/Conjugate_prior#When_likelihood_function_is_a_continuous_distribution

[15] симметрии: http://www.braintools.ru/article/3088

[16] ScipyTTestInd: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html

[17] Chi2Test: https://en.wikipedia.org/wiki/Chi-squared_test

[18] Chi2Pearson: https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test

[19] Chi2Dist: https://en.wikipedia.org/wiki/Chi-squared_distribution

[20] ScipyChi2Con: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2_contingency.html

[21] MannWhitneyU: https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test

[22] NormalSum: https://en.wikipedia.org/wiki/Sum_of_normally_distributed_random_variables

[23] ScipyMannWhitneyU: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.mannwhitneyu.html

[24] Источник: https://habr.com/ru/articles/987764/?utm_campaign=987764&utm_source=habrahabr&utm_medium=rss

www.BrainTools.ru

Rambler's Top100