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

Байесовский подход применен к А/Б-тесту конверсий с 3 группами. Лучшая группа выбирается сравнением апостериорных распределений. Способ применим для других метрик и большего количества вариантов.

Блокнот: https://github.com/andrewbrdk/Bayesian-AB-Testing/blob/main/appendices/Множественные_сравнения.ipynb .

Библиотеки

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

np.random.seed(7)

В А/Б-тестах бывает больше 2 вариантов. "Проверка статистических гипотез" в таких случаях требует поправок на множественные сравнения [MultipleComp, FWER, Bonf]. В байесовском подходе лучшая группа выбирается сравнением апостериорных распределений, дополнительные поправки не требуются.

На три версии веб-страницы A, B и С зашло по N=1000 человек. Кнопку "Продолжить" нажали $n_{s_A}=100$ , $n_{s_B}=105$ , $n_{s_C}=110$ человек соответственно. С какой вероятностью конверсия каждого из вариантов лучшая?

Для каждого варианта нужно оценить вероятность наибольшей конверсии из всех групп: $P(\text{Best } A) \equiv P(p_A > p_B \cap p_A > p_C)$ для А, аналогично для B и C. Вероятности можно оценить численно сравнением выборок апостериорных распределений. Для конверсий правдоподобие $P(\mathcal{D} | \mathcal{H})$ задается биномиальным распределением, априорное распределение $P(\mathcal{H})$ - бета-распределением. В таком случае апостериорные распределения $P(\mathcal{H} | \mathcal{D})$ также будут бета-распределениями [BayesABConv, BetaDist, SciPyBeta, ConjPrior].

$P(\mathcal{H} | \mathcal{D}) \propto P(\mathcal{D} | \mathcal{H}) P(\mathcal{H})$ $P(\mathcal{D} | \mathcal{H}) = P(n_s, N | p) = \mbox{Binom}(n_s, N | p) = C_{N}^{n_s} p^{n_s} (1-p)^{N-n_s}$ $P(\mathcal{H}) = P(p) = \mbox{Beta}(p; \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} p^{\alpha-1}(1-p)^{\beta-1}$ $\begin{split}P(\mathcal{H} | \mathcal{D}) & = P(p | n_s, N) = \mbox{Beta}(p; \alpha + n_s, \beta + N - n_s)\end{split}$

На графике приведены апостериорные распределения каждой группы. Распределения пересекаются. Вероятности лучшей конверсии $P(\text{Best } A) = 15\%$ , $P(\text{Best } B) = 30\%$ , $P(\text{Best } C) = 55\%$ .

Апостериорные распределения конверсий

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)

N = 1000
sa = 100
sb = 105
sc = 110

p_dist_a = posterior_dist_binom(ns=sa, ntotal=N)
p_dist_b = posterior_dist_binom(ns=sb, ntotal=N)
p_dist_c = posterior_dist_binom(ns=sc, ntotal=N)

npost = 50000
samp_a = p_dist_a.rvs(size=npost)
samp_b = p_dist_b.rvs(size=npost)
samp_c = p_dist_c.rvs(size=npost)

p_a_best = np.sum((samp_a > samp_b) & (samp_a > samp_c)) / npost
p_b_best = np.sum((samp_b > samp_a) & (samp_b > samp_c)) / npost
p_c_best = np.sum((samp_c > samp_a) & (samp_c > samp_b)) / npost

xaxis_max = 0.2
x = np.linspace(0, xaxis_max, 1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=p_dist_a.pdf(x), line_color='black', name='A'))
fig.add_trace(go.Scatter(x=x, y=p_dist_b.pdf(x), line_color='black', line_dash='longdash', name='B'))
fig.add_trace(go.Scatter(x=x, y=p_dist_c.pdf(x), line_color='black', line_dash='dot', name='C'))
fig.update_layout(title='Апостериорные распределения',
                  xaxis_title='$p$',
                  yaxis_title='Плотность вероятности',
                  xaxis_range=[0, xaxis_max],
                  hovermode="x",
                  height=500)
fig.show()

print(f"P Best:")
print(f"P(Best A) = P(A>B & A>C) = {p_a_best}")
print(f"P(Best B) = P(B>A & B>C) = {p_b_best}")
print(f"P(Best C) = P(C>A & C>B) = {p_c_best}")

P Best:
P(Best A) = P(A>B & A>C) = 0.14664
P(Best B) = P(B>A & B>C) = 0.30228
P(Best C) = P(C>A & C>B) = 0.55108

Апостериорные распределения конверсий 3 групп. Вероятности лучшей конверсии P(Best A) = 15%, P(Best B) = 30%, P(Best C) = 55%.

Количество правильно угаданных вариантов в серии экспериментов следующее. В группе A задается конверсия p = 0.1, в группах B и C конверсия выбирается случайно в диапазоне $\pm 5\%$ от p. В группах генерируются данные с шагом n_samp_step. На каждом шаге в каждом варианте считаются апостериорные распределения и вероятность лучшей конверсии среди всех групп $P(\text{Best } A)$ и др. Эксперимент останавливается, если в одной из групп вероятность наибольшей конверсии достигает prob_stop=0.95 или сгенерировано максимальное количество точек n_samp_max. Проводится nexps экспериментов, считается доля правильно угаданных групп. В данном случае в nexps = 1000 правильно угадано 951. Точность 0.951 близка ожидаемой prob_stop = 0.95.

Nexp: 1000, Correct Guesses: 951, Accuracy: 0.951

Правильно угаданные варианты в серии экспериментов

def p_best(*args, n_post_samp=50_000):
    samp = [d.rvs(size=n_post_samp) for d in args]
    best_group = np.argmax(np.vstack(samp), axis=0)
    u = np.unique(best_group, return_counts=True)
    p_best = np.zeros(len(args))
    for i, c in zip(u[0], u[1]):
        p_best[i] = c
    p_best = p_best / n_post_samp
    return p_best

cmp = pd.DataFrame(columns=['A', 'B', 'C', 'best_exact', 
                            'exp_samp_size', 'A_exp', 'B_exp', 'C_exp', 
                            'best_exp', 'p_best'])

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['C'] = p * (1 + stats.uniform.rvs(loc=-0.05, scale=0.1, size=nexps))
cmp['best_exact'] = cmp.apply(lambda r: 'A' if r['A'] > r['B'] and r['A'] > r['C'] else 'B' if r['B'] > r['A'] and r['B'] > r['C'] else 'C', axis=1)

n_samp_max = 30_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']
    pC = cmp.at[i, 'C']
    exact_dist_A = stats.bernoulli(p=pA)
    exact_dist_B = stats.bernoulli(p=pB)
    exact_dist_C = stats.bernoulli(p=pC)
    n_samp_total = 0
    ns_A = 0
    ns_B = 0
    ns_C = 0
    while n_samp_total < n_samp_max:
        dA = exact_dist_A.rvs(n_samp_step)
        dB = exact_dist_B.rvs(n_samp_step)
        dC = exact_dist_C.rvs(n_samp_step)
        n_samp_total += n_samp_step
        ns_A = ns_A + np.sum(dA)
        ns_B = ns_B + np.sum(dB)
        ns_C = ns_C + np.sum(dC)
        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)
        post_dist_C = posterior_dist_binom(ns=ns_C, ntotal=n_samp_total)
        p_best_A, p_best_B, p_best_C = p_best(post_dist_A, post_dist_B, post_dist_C)
        best_gr = 'A' if p_best_A >= prob_stop else 'B' if  p_best_B >= prob_stop else 'C' if p_best_C >= prob_stop else None
        if best_gr:
            cmp.at[i, 'A_exp'] = post_dist_A.mean()
            cmp.at[i, 'B_exp'] = post_dist_B.mean()
            cmp.at[i, 'C_exp'] = post_dist_C.mean()
            cmp.at[i, 'exp_samp_size'] = n_samp_total
            cmp.at[i, 'best_exp'] = best_gr
            cmp.at[i, 'p_best'] = max(p_best_A, p_best_B, p_best_C)
            break
    print(f'done {i}: nsamp {n_samp_total}, best_gr {best_gr}, P_best {max(p_best_A, p_best_B, p_best_C)}')

cmp['correct'] = cmp['best_exact'] == cmp['best_exp']
display(cmp.head(20))
cor_guess = np.sum(cmp['correct'])
print(f"Nexp: {nexps}, Correct Guesses: {cor_guess}, Accuracy: {cor_guess / nexps}")

Ссылки:

[BayesABConv] - Bayesian A/B-Testing, GitHub.
[BetaDist] - Beta Distribution, Wikipedia.
[Bonf] - Bonferroni Correction, Wikipedia.
[ConjPrior] - Conjugate Prior, Wikipedia.
[FWER] - Family-wise Error Rate, Wikipedia.
[MultipleComp] - Multiple Comparisons Problem, Wikipedia.
[SciPyBeta] - scipy.stats.beta, SciPy Reference.