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

import numpy as np import pandas as pd import scipy.stats as stats import plotly.graph_objects as go  np.random.seed(7)

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

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