A/Bテストの統計学 — 仮説検定とベイズ的アプローチ

import numpy as np
from scipy import stats

def ab_test_z(x_a, n_a, x_b, n_b, alpha=0.05):
    """比率の差のZ検定（両側）"""
    # 標本比率
    p_a = x_a / n_a
    p_b = x_b / n_b

    # プール推定量
    p_pool = (x_a + x_b) / (n_a + n_b)

    # 検定統計量
    se = np.sqrt(p_pool * (1 - p_pool) * (1/n_a + 1/n_b))
    z = (p_b - p_a) / se

    # p値
    p_value = 2 * stats.norm.sf(abs(z))

    # 信頼区間（非プール分散を使用）
    se_ci = np.sqrt(p_a*(1-p_a)/n_a + p_b*(1-p_b)/n_b)
    z_crit = stats.norm.ppf(1 - alpha/2)
    ci_lower = (p_b - p_a) - z_crit * se_ci
    ci_upper = (p_b - p_a) + z_crit * se_ci

    return {
        'p_a': p_a, 'p_b': p_b,
        'diff': p_b - p_a,
        'z_stat': z, 'p_value': p_value,
        'ci': (ci_lower, ci_upper),
        'significant': p_value < alpha
    }

# 例: A群 1000人中100人コンバージョン、B群 1000人中130人コンバージョン
result = ab_test_z(100, 1000, 130, 1000)
print("=== 頻度主義的 A/Bテスト ===")
print(f"コンバージョン率: A = {result['p_a']:.3f}, B = {result['p_b']:.3f}")
print(f"差: {result['diff']:.3f}")
print(f"Z統計量: {result['z_stat']:.3f}")
print(f"p値: {result['p_value']:.4f}")
print(f"95%信頼区間: ({result['ci'][0]:.4f}, {result['ci'][1]:.4f})")
print(f"統計的に有意: {result['significant']}")

def sample_size_two_proportions(p_a, mde, alpha=0.05, power=0.80, sides=2):
    """2群の比率の差のA/Bテストに必要なサンプルサイズを計算"""
    p_b = p_a + mde

    if sides == 2:
        z_alpha = stats.norm.ppf(1 - alpha/2)
    else:
        z_alpha = stats.norm.ppf(1 - alpha)
    z_beta = stats.norm.ppf(power)

    # 各群のサンプルサイズ
    n = ((z_alpha + z_beta)**2 * (p_a*(1-p_a) + p_b*(1-p_b))) / mde**2
    return int(np.ceil(n))

# 様々な条件でのサンプルサイズ
print("必要サンプルサイズ（各群、α=0.05, Power=0.80, 両側）:")
print("-" * 50)
for p_a in [0.05, 0.10, 0.20]:
    for rel_mde in [0.05, 0.10, 0.20]:
        mde = p_a * rel_mde
        n = sample_size_two_proportions(p_a, mde)
        print(f"  p_A={p_a:.2f}, 相対MDE={rel_mde:.0%} (絶対MDE={mde:.4f}): "
              f"n={n:,} /群, 合計={2*n:,}")

def bayesian_ab_test(x_a, n_a, x_b, n_b,
                     prior_alpha=1, prior_beta=1,
                     n_samples=100000):
    """ベイズ的A/Bテスト"""
    # 事後分布のパラメータ
    post_a_alpha = prior_alpha + x_a
    post_a_beta = prior_beta + n_a - x_a
    post_b_alpha = prior_alpha + x_b
    post_b_beta = prior_beta + n_b - x_b

    # 事後分布からサンプリング
    samples_a = np.random.beta(post_a_alpha, post_a_beta, n_samples)
    samples_b = np.random.beta(post_b_alpha, post_b_beta, n_samples)

    # Bが勝つ確率
    prob_b_wins = np.mean(samples_b > samples_a)

    # 差の分布
    diff_samples = samples_b - samples_a

    # 期待損失
    loss_b = np.mean(np.maximum(samples_a - samples_b, 0))  # Bを選んだときの損失
    loss_a = np.mean(np.maximum(samples_b - samples_a, 0))  # Aを選んだときの損失

    # 相対的な改善率の分布
    relative_uplift = (samples_b - samples_a) / samples_a

    return {
        'prob_b_wins': prob_b_wins,
        'diff_mean': np.mean(diff_samples),
        'diff_ci': (np.percentile(diff_samples, 2.5),
                    np.percentile(diff_samples, 97.5)),
        'loss_a': loss_a,
        'loss_b': loss_b,
        'relative_uplift_mean': np.mean(relative_uplift),
        'relative_uplift_ci': (np.percentile(relative_uplift, 2.5),
                               np.percentile(relative_uplift, 97.5)),
        'samples_a': samples_a,
        'samples_b': samples_b,
        'diff_samples': diff_samples
    }

np.random.seed(42)
bayes = bayesian_ab_test(100, 1000, 130, 1000)
print("\n=== ベイズ的 A/Bテスト ===")
print(f"P(B > A | data) = {bayes['prob_b_wins']:.4f}")
print(f"差の事後平均: {bayes['diff_mean']:.4f}")
print(f"差の95%信用区間: ({bayes['diff_ci'][0]:.4f}, {bayes['diff_ci'][1]:.4f})")
print(f"期待損失 (Bを選択): {bayes['loss_b']:.5f}")
print(f"期待損失 (Aを選択): {bayes['loss_a']:.5f}")
print(f"相対的改善率の事後平均: {bayes['relative_uplift_mean']:.2%}")
print(f"相対的改善率の95%信用区間: ({bayes['relative_uplift_ci'][0]:.2%}, "
      f"{bayes['relative_uplift_ci'][1]:.2%})")

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# テストデータ
x_a, n_a = 100, 1000
x_b, n_b = 130, 1000

fig, axes = plt.subplots(2, 3, figsize=(18, 11))

# (a) 事後分布の比較
ax = axes[0][0]
post_a = stats.beta(1 + x_a, 1 + n_a - x_a)
post_b = stats.beta(1 + x_b, 1 + n_b - x_b)
x_range = np.linspace(0.06, 0.18, 500)

ax.plot(x_range, post_a.pdf(x_range), 'b-', linewidth=2, label=f'A posterior (CVR={x_a/n_a:.1%})')
ax.fill_between(x_range, post_a.pdf(x_range), alpha=0.2, color='blue')
ax.plot(x_range, post_b.pdf(x_range), 'r-', linewidth=2, label=f'B posterior (CVR={x_b/n_b:.1%})')
ax.fill_between(x_range, post_b.pdf(x_range), alpha=0.2, color='red')

ax.set_xlabel('Conversion rate', fontsize=11)
ax.set_ylabel('Density', fontsize=11)
ax.set_title('(a) Posterior Distributions', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) 差の事後分布
ax = axes[0][1]
samples_a = np.random.beta(1 + x_a, 1 + n_a - x_a, 100000)
samples_b = np.random.beta(1 + x_b, 1 + n_b - x_b, 100000)
diff = samples_b - samples_a

ax.hist(diff, bins=100, density=True, alpha=0.7, color='green', edgecolor='none')
ax.axvline(0, color='black', linewidth=2, linestyle='--', label='No difference')
ax.axvline(np.mean(diff), color='red', linewidth=2, label=f'Mean = {np.mean(diff):.4f}')

ci_low, ci_high = np.percentile(diff, [2.5, 97.5])
ax.axvline(ci_low, color='orange', linewidth=1.5, linestyle=':', label=f'95% CI: [{ci_low:.4f}, {ci_high:.4f}]')
ax.axvline(ci_high, color='orange', linewidth=1.5, linestyle=':')

prob_positive = np.mean(diff > 0)
ax.fill_between([0, max(diff)], 0, 0.1, alpha=0.1, color='green',
                transform=ax.get_xaxis_transform())
ax.set_xlabel('$p_B - p_A$', fontsize=11)
ax.set_ylabel('Density', fontsize=11)
ax.set_title(f'(b) Posterior of Difference (P(B>A)={prob_positive:.1%})', fontsize=12)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# (c) 覗き見問題のシミュレーション
ax = axes[0][2]
n_sims = 2000
max_n = 5000
check_points = np.arange(100, max_n + 1, 50)
true_p = 0.10  # 帰無仮説が真（両群とも同じ比率）

false_positive_fixed = []
false_positive_peeking = []

for _ in range(n_sims):
    data_a = np.random.binomial(1, true_p, max_n)
    data_b = np.random.binomial(1, true_p, max_n)

    peeking_rejected = False
    for cp in check_points:
        xa = data_a[:cp].sum()
        xb = data_b[:cp].sum()
        p_pool = (xa + xb) / (2 * cp)
        if p_pool == 0 or p_pool == 1:
            continue
        se = np.sqrt(p_pool * (1-p_pool) * 2/cp)
        z = (xb/cp - xa/cp) / se
        p_val = 2 * stats.norm.sf(abs(z))
        if p_val < 0.05:
            peeking_rejected = True
            break

    false_positive_peeking.append(peeking_rejected)

    # 固定サンプルサイズでの検定
    xa = data_a.sum()
    xb = data_b.sum()
    p_pool = (xa + xb) / (2 * max_n)
    se = np.sqrt(p_pool * (1-p_pool) * 2/max_n)
    z = (xb/max_n - xa/max_n) / se
    p_val = 2 * stats.norm.sf(abs(z))
    false_positive_fixed.append(p_val < 0.05)

fpr_peeking = np.mean(false_positive_peeking)
fpr_fixed = np.mean(false_positive_fixed)

# 確認回数による偽陽性率の変化
check_counts = [1, 2, 5, 10, 20, 50, 100]
fprs = []
for k in check_counts:
    cp_subset = np.linspace(100, max_n, k, dtype=int)
    fp_count = 0
    for _ in range(n_sims):
        da = np.random.binomial(1, true_p, max_n)
        db = np.random.binomial(1, true_p, max_n)
        for cp in cp_subset:
            xa = da[:cp].sum()
            xb = db[:cp].sum()
            pp = (xa + xb) / (2 * cp)
            if pp == 0 or pp == 1:
                continue
            se = np.sqrt(pp * (1-pp) * 2/cp)
            z = (xb/cp - xa/cp) / se
            if 2 * stats.norm.sf(abs(z)) < 0.05:
                fp_count += 1
                break
    fprs.append(fp_count / n_sims)

ax.plot(check_counts, fprs, 'ro-', linewidth=2, markersize=6, label='Peeking FPR')
ax.axhline(0.05, color='blue', linewidth=2, linestyle='--', label='Nominal $\\alpha = 0.05$')
ax.set_xlabel('Number of peeks', fontsize=11)
ax.set_ylabel('False Positive Rate', fontsize=11)
ax.set_title('(c) Peeking Problem', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xscale('log')

# (d) サンプルサイズと必要期間
ax = axes[1][0]
daily_traffic = [100, 500, 1000, 5000, 10000, 50000]
p_a_base = 0.10
rel_mdes = [0.05, 0.10, 0.20]
z_a = stats.norm.ppf(0.975)
z_b = stats.norm.ppf(0.80)

bar_width = 0.25
x_pos = np.arange(len(daily_traffic))

for i, rel_mde in enumerate(rel_mdes):
    mde = p_a_base * rel_mde
    p_b = p_a_base + mde
    n_per_group = ((z_a + z_b)**2 * (p_a_base*(1-p_a_base) + p_b*(1-p_b))) / mde**2
    n_per_group = int(np.ceil(n_per_group))

    days_needed = [2 * n_per_group / dt for dt in daily_traffic]
    ax.bar(x_pos + i*bar_width, days_needed, bar_width,
           label=f'Relative MDE = {rel_mde:.0%}', alpha=0.8)

ax.set_xlabel('Daily traffic per group', fontsize=11)
ax.set_ylabel('Days needed', fontsize=11)
ax.set_title(f'(d) Test Duration ($p_A$={p_a_base:.0%})', fontsize=12)
ax.set_xticks(x_pos + bar_width)
ax.set_xticklabels([f'{d:,}' for d in daily_traffic], rotation=30)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis='y')
ax.set_yscale('log')

# (e) ベイズ的A/Bテストの経時変化
ax = axes[1][1]
true_pa = 0.10
true_pb = 0.13
max_obs = 2000
np.random.seed(123)

data_a_seq = np.random.binomial(1, true_pa, max_obs)
data_b_seq = np.random.binomial(1, true_pb, max_obs)

obs_points = np.arange(20, max_obs + 1, 10)
prob_b_wins_seq = []
loss_b_seq = []

for t in obs_points:
    xa = data_a_seq[:t].sum()
    xb = data_b_seq[:t].sum()
    sa = np.random.beta(1 + xa, 1 + t - xa, 50000)
    sb = np.random.beta(1 + xb, 1 + t - xb, 50000)
    prob_b_wins_seq.append(np.mean(sb > sa))
    loss_b_seq.append(np.mean(np.maximum(sa - sb, 0)))

ax.plot(obs_points, prob_b_wins_seq, 'r-', linewidth=2, label='$P(p_B > p_A | \\mathrm{data})$')
ax.axhline(0.95, color='green', linewidth=1.5, linestyle='--', label='Threshold = 0.95')
ax.axhline(0.5, color='gray', linewidth=1, linestyle=':', alpha=0.5)

# 95%を初めて超えた時点
first_cross = None
for i, p in enumerate(prob_b_wins_seq):
    if p >= 0.95:
        first_cross = obs_points[i]
        break
if first_cross:
    ax.axvline(first_cross, color='green', linewidth=1, linestyle='-', alpha=0.3)
    ax.annotate(f'First crosses 95%\nat n={first_cross}',
                xy=(first_cross, 0.95), fontsize=9,
                xytext=(first_cross + 200, 0.75),
                arrowprops=dict(arrowstyle='->', color='green'))

ax.set_xlabel('Observations per group', fontsize=11)
ax.set_ylabel('Probability', fontsize=11)
ax.set_title(f'(e) Bayesian A/B Test Over Time\n($p_A$={true_pa}, $p_B$={true_pb})', fontsize=12)
ax.legend(fontsize=9, loc='center right')
ax.grid(True, alpha=0.3)
ax.set_ylim(0.3, 1.05)

# (f) 頻度主義 vs ベイズの判定結果比較
ax = axes[1][2]
np.random.seed(42)
n_experiments = 500
p_a_true = 0.10
effects = np.random.choice([0, 0.01, 0.02, 0.03, 0.05], n_experiments,
                           p=[0.4, 0.15, 0.15, 0.15, 0.15])
test_n = 500

freq_decisions = []
bayes_decisions = []

for eff in effects:
    p_b_true = p_a_true + eff
    xa = np.random.binomial(test_n, p_a_true)
    xb = np.random.binomial(test_n, p_b_true)

    # 頻度主義
    p_pool = (xa + xb) / (2 * test_n)
    if p_pool == 0 or p_pool == 1:
        freq_decisions.append(False)
    else:
        se = np.sqrt(p_pool * (1-p_pool) * 2 / test_n)
        z = (xb/test_n - xa/test_n) / se
        freq_decisions.append(2 * stats.norm.sf(abs(z)) < 0.05)

    # ベイズ（期待損失 < 0.001 かつ P(B>A) > 0.95）
    sa = np.random.beta(1 + xa, 1 + test_n - xa, 50000)
    sb = np.random.beta(1 + xb, 1 + test_n - xb, 50000)
    p_b_better = np.mean(sb > sa)
    exp_loss = np.mean(np.maximum(sa - sb, 0))
    bayes_decisions.append(p_b_better > 0.95 and exp_loss < 0.001)

# 効果あり/なし別の判定率
true_null = effects == 0
true_alt = effects > 0

results = {
    'Frequentist': {
        'FPR': np.mean(np.array(freq_decisions)[true_null]),
        'TPR': np.mean(np.array(freq_decisions)[true_alt]),
    },
    'Bayesian': {
        'FPR': np.mean(np.array(bayes_decisions)[true_null]),
        'TPR': np.mean(np.array(bayes_decisions)[true_alt]),
    }
}

categories = ['FPR\n(False Positive)', 'TPR\n(True Positive)']
freq_vals = [results['Frequentist']['FPR'], results['Frequentist']['TPR']]
bayes_vals = [results['Bayesian']['FPR'], results['Bayesian']['TPR']]

x_cat = np.arange(len(categories))
ax.bar(x_cat - 0.15, freq_vals, 0.3, label='Frequentist', color='steelblue', alpha=0.8)
ax.bar(x_cat + 0.15, bayes_vals, 0.3, label='Bayesian', color='coral', alpha=0.8)

for i, (fv, bv) in enumerate(zip(freq_vals, bayes_vals)):
    ax.text(i - 0.15, fv + 0.01, f'{fv:.2f}', ha='center', fontsize=10)
    ax.text(i + 0.15, bv + 0.01, f'{bv:.2f}', ha='center', fontsize=10)

ax.set_xticks(x_cat)
ax.set_xticklabels(categories, fontsize=11)
ax.set_ylabel('Rate', fontsize=11)
ax.set_title(f'(f) Frequentist vs Bayesian (n={test_n}/group)', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis='y')
ax.set_ylim(0, max(max(freq_vals), max(bayes_vals)) + 0.15)

plt.tight_layout()
plt.savefig('ab_testing_statistics.png', dpi=150, bbox_inches='tight')
plt.show()

def check_srm(n_a, n_b, expected_ratio=0.5):
    """サンプル比率の不一致を検出"""
    n_total = n_a + n_b
    expected_a = n_total * expected_ratio
    expected_b = n_total * (1 - expected_ratio)

    chi2 = (n_a - expected_a)**2 / expected_a + (n_b - expected_b)**2 / expected_b
    p_value = 1 - stats.chi2.cdf(chi2, df=1)

    return {
        'chi2': chi2,
        'p_value': p_value,
        'ratio': n_a / n_total,
        'srm_detected': p_value < 0.001  # 厳しめの閾値
    }

# 正常な場合
srm1 = check_srm(4980, 5020)
print(f"\n=== SRM検出 ===")
print(f"正常例 (4980 vs 5020): p={srm1['p_value']:.4f}, SRM={srm1['srm_detected']}")

# 異常な場合
srm2 = check_srm(4700, 5300)
print(f"異常例 (4700 vs 5300): p={srm2['p_value']:.6f}, SRM={srm2['srm_detected']}")