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

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

- P-значения
- Т-тест
- Тест $\chi^2$
- U-критерий Манна-Уитни
- Ссылки

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

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

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

Решение о гипотезе принимают по -значению $p = P_{T}(x \ge x_{0} | H_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}$

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

Т-тест

Средние сравнивают -тестами [TTest]. Пусть есть выборки размера из двух случайных величин . В предположении одинаковых точных средних для отношения разности выборочных средних к стандартной ошибке разности средних $X = \overline{\Delta}/s_{\Delta}$ ожидают -распределение [WelchT]. При достаточно большом количестве данных оно близко стандартному нормальному $\text{Norm}(0, 1)$ [TDist]. По выборкам считают фактическое отношение $x_0 = \overline{\Delta}/s_{\Delta}$ . Вычисляют вероятность получить или большее значение - одностороннее -значение $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_{X}(x \ge x_0 | H_0)$ интересна вероятность среднего больше при условии собранных данных $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]. При достаточно широком априорном распределении $\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}$

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

$\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)$ сравнивается с -значением -теста. Используется односторонний -тест с разными дисперсиями групп (equal_var=False, alternative) [ScipyTTestInd]. Вероятность
$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)$ закрашена светлым. По свойствам нормального распределения площади этих областей совпадают. Таким образом -значение численно близко байесовской оценке вероятности. Стоит помнить, что они не эквивалентны - у них разные определения.

Т-распределение и апостериорное среднего разности

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'$T\text{-распределение и апостериорное среднего разности}$',
                  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(\mu_{\Delta} > 0 | \overline{\Delta})$ можно проверить по количеству правильно угаданных вариантов в серии экспериментов. В каждом эксперименте задается два нормальных распределения. В группе среднее фиксировано mu=0.1, в выбирается случайно в диапазоне $\pm 5\%$ от mu. В группы добавляются данные по n_samp_step точек за шаг. На каждом шаге считается -тест. Эксперимент останавливается, если достигает prob_stop=0.95 или использовано максимальное количество точек n_samp_max. Длительность эксперимента не фиксируется заранее. Миниальный размер выборки n_samp_min + n_samp_step. При остановке эксперимента для сравнения с -значением считается байесовское апострериорное распределение и вероятность $P(\mu_{\Delta} > 0 | \overline{\Delta})$ . Процедура повторяется nexps раз, считается доля правильно угаданных групп во всех экспериментах. Из nexps = 1000 завершено 880, в них правильно угадано 818 вариантов. Точность 0.93 близка prob_stop = 0.95. -значения в каждом эксперименте близки байесовским вероятностям. При недостаточном 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}")

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

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

Для А/Б-теста конверсий с двумя группами в предположении одинаковых ожидаемых конверсий тестовую статистику можно привести к виду $\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_{\chi_1^2}(x > \chi^2 | H_0)$ соответствует областям $P\left(\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 P\left( \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)$ . При типичных значениях бета-распределение близко нормальному. Распределение разности конверсий также будет нормальным. При равномерных априорных распределениях и слабой разнице между группами $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 P\left(\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^2\right) \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 P\left(\text{Norm}(x > 0; \Delta p, s_{\Delta}^2) \right) = P\left(\text{Norm}(x > 0; \chi, 1) \right), \quad \chi = \Delta p / s_{\Delta} \end{split}$

По симметрии $P\left( \text{Norm}(x > \chi | 0, 1) \right) = P\left(\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 P\left( \text{Norm}(x > \chi; 0, 1) \right) \\ & = 2 P\left(\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$ проверяется по выборке из двух распределений Бернулли с конверсиями и . Данные для $\chi^2$ -теста задаются в виде таблицы со строками и [ScipyChi2Con]. -значение $\chi^2$ -теста не позволяет выбрать между и , поэтому дополнительно сравниваются средние в выборках. Распределение $\chi^2_1$ на первом графике ниже. Закрашенная область соответствует -значению $\mbox{p-val} = P_{\chi_1^2}(x > \chi^2 | H_0)$ . На втором графике закрашенные темные области $x > \chi$ и $x < - \chi$ соответствуют области -значения при возведении в квадрат. Серый график - нормальное распределение $\text{Norm}(\chi, 1)$ . Закрашенная область серого графика приближенно равна $1 - P(\theta_B > \theta_A | A_i, B_j)$ . Площади закрашенной серой и каждой из темных областей совпадают. Связь -значения и байесовской оценки вероятности численно выполняется.

Распределение Хи-квадрат, 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'$P\text{-значение и байесовская модель}$',
                  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.

-значение $\chi^2$ -теста используется для выбора вариантов с большей конверсией в серии экспериментов. В каждом эксперименте 2 группы, конверсия фиксирована, выбирается случайно в диапазоне $\pm5\%$ от . В каждой группе добавляются данные по n_samp_step точек за шаг. На каждом шаге вычисляется -значение $\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}")

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

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

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

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

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

В байесовском подходе вероятность можно оценить сравнением апостериорных предиктивных распределений. Для этого нужно задать модели исходных распределений, построить апостериорные распределения параметров, апостериорные предиктивные распределения и по ним оценить . Другой вариант, не требующий предположений о распределениях - моделировать вероятность $\theta = P(B > A)$ точки больше . Данными будут пары , в каждой паре точка больше с вероятностью $\theta$ , по общему количеству пар и количеству пар с больше нужно оценить $\theta$ . Такая задача похожа на моделирование конверсий. Если при формировании пар из исходных выборок каждую точку сравнивать с каждой точкой , пары с одинаковыми или будут зависимы. Для независимости каждая точка сравнивается только с одной точкой . Аналогично конверсиям правдоподобие задается биномиальным распределением $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}$

В -статистике каждая точка сравнивается с каждой , в байесовской модели точки участвуют только в одном сравнении. Из-за меньшего использования данных дисперсия байесовской модели шире -статистики. -значение не будет близко совпадать с байесовской вероятностью $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}$

Проверка связи -значения с байесовской вероятностью сделана на примере нормальных распределений. Вероятность точки из распределения больше равна вероятности разности больше нуля . Разность нормальных распределений также нормальное распределение со средним, равным разности исходных средних, и дисперсией, равной сумме дисперсий [NormalSum].

$\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 и единичной дисперсией. Темная линия на втором графике - в предположении эквивалентности распределений. Темная закрашенная область соответствует -значению [ScipyMannWhitneyU]. Светлая линия показывает байесовскую модель, светлая закрашенная область - байесовскую вероятность . Байесовская модель шире из-за большей дисперсии. Площадь светлой области обычно больше темной, что соответствует $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'$U\text{-статистика и байесовская модель}$',
    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.

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

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}")

Для проверки правильно угаданных вариантов по -статистике выставлен больший шаг - добавляется по 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}")

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

Ссылки

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