贝叶斯检验

🧠 一、什么是贝叶斯检验

用一句话概括：贝叶斯检验就是"带着常识（先验），看了证据（数据）后，更新你的信念（后验），并用一个叫贝叶斯因子（BF）的刻度尺衡量证据有多硬"。

1. 先验概率（Prior）

在看到数据之前，你对世界的"朴素想法"或常识。例如，在进行 A/B 测试之前，我们通常不会认为 B 版本一定远远优于 A 版本。

可能认为：

• A 和 B 效果差不多；

• B 可能略好一点；

• B 可能甚至更差。

这些已有经验，就是先验。

2. 似然（Likelihood）

在某个假设下，当前数据出现的可能性有多大。

例如假设：

• H0：A 版本和 B 版本没有区别；

• H1：B 版本效果更好（或更一般地：A、B 各自独立，不强制相等）。

那么我们需要比较：如果 H0 成立，当前实验结果是否合理？如果 H1 成立，当前实验结果是否合理？

3. 后验概率（Posterior）

看到数据之后，你更新后的新理解。

贝叶斯公式：

其中 P(H) 是先验概率，P(D|H) 是似然，P(H|D) 是后验概率。

4. 贝叶斯因子（Bayes Factor, BF）

衡量"数据更支持哪个假设"的证据强度，贝叶斯检验中最重要的指标：

它表示当前数据支持 H1 的程度，是支持 H0 的多少倍。经验参考（Jeffreys 分级）：

• BF ≈ 1：数据无法区分两个假设；

• BF > 3：开始支持 H1；

• BF > 10：强烈支持 H1；

• 反过来，BF < 1（比如 0.1）就意味着数据反而在支持 H0，而不是"没有证据"。这一点很容易被忽略，但在后面的实战部分非常关键。

⚖️ 二、贝叶斯检验 vs 传统 P 值检验

🍬 三、通俗案例-糖果罐

假设你有两罐糖：

• 罐子 A：红糖 50% / 蓝糖 50%

• 罐子 B：红糖 80% / 蓝糖 20%

场景：你蒙眼抓了 5 颗糖，结果是 4 红 1 蓝。问：这更可能是哪个罐子？

计算过程：

1. 如果是罐子 A：抓到"4 红 1 蓝"的概率 ≈ 0.156

2. 如果是罐子 B：抓到"4 红 1 蓝"的概率 ≈ 0.410

结论：证据比（贝叶斯因子近似）为 0.410/0.156 ≈ 2.6。数据对"这是罐子 B"的支持度是"罐子 A"的 2.6 倍。虽然不是绝对铁证，但明显更偏向 B。

🛠️ 四、A/B 测试全流程

1. 核心数学工具

在 A/B 转化率（伯努利–贝塔共轭）场景下，可以用闭式解快速计算，无需 MCMC 采样。

• 贝叶斯因子公式：（B 为 Beta 函数）

• ROPE（实用等价区间）：定义一个业务上无感的差异范围，看后验落在该区间内的概率，判断是否"业务等价"。这个区间不是随手拍的——比如可以按"多大的转化率差异，换算成的收入增量不足以覆盖切换/维护成本"来倒推，本文取 ±0.2 个百分点仅作演示。

2. 先验的选择

一个很容易踩的坑：如果直接用 Beta(2,2) 之类"看起来温和"的先验，它的均值和众数其实都在 0.5。但电商/APP 场景的转化率通常只有百分之几，这种先验相当于先入为主地认为"转化率大概率落在 30%~70%附近"，和常识明显不符，只是样本量够大时会被数据"压过去"，不容易被发现。

更合理的做法是：围绕业务历史基线转化率展开先验，再用一个"强度"参数控制先验有多固执：

def beta_params_from_mean_strength(mean, strength):
    """按'先验均值 + 先验强度(=a+b)'两个可解释的业务维度反推 Beta 参数"""
    a = mean * strength
    b = (1 - mean) * strength
    return a, b

# 假设历史基线转化率约 5%，先验强度取 4（约等于"4次历史观测"的信息量）
PRIOR_MEAN = 0.05
PRIOR_STRENGTH = 4
a_prior, b_prior = beta_params_from_mean_strength(PRIOR_MEAN, PRIOR_STRENGTH)

3. 生成虚拟数据

模拟 A/B 两组转化数据（A 组真值 5%，B 组真值 6%）。

import numpy as np
import pandas as pd
from scipy.stats import beta as beta_dist
from scipy.special import betaln

np.random.seed(42)
p_true_A, p_true_B = 0.050, 0.060
n_A, n_B = 5000, 4800

A = np.random.binomial(1, p_true_A, size=n_A)
B = np.random.binomial(1, p_true_B, size=n_B)
x_A, x_B = A.sum(), B.sum()

print(f"观测: A 组 {x_A}/{n_A}, B 组 {x_B}/{n_B}")

4. 计算贝叶斯因子

def bayes_factor_beta_binomial(xA, nA, xB, nB, a=1, b=1):
    # 默认 Beta(1,1) 无信息先验；正式分析请显式传入业务先验，
    # 避免"默认值"和"正文实际使用值"不一致造成误解
    num = betaln(xA + a, nA - xA + b) + betaln(xB + a, nB - xB + b)
    den = betaln(a, b) + betaln(xA + xB + a, (nA + nB) - (xA + xB) + b)
    return np.exp(num - den)

BF10 = bayes_factor_beta_binomial(x_A, n_A, x_B, n_B, a=a_prior, b=b_prior)
print(f"BF10 = {BF10:.4f}  (BF01 = {1/BF10:.2f})")

跑出来的结果是：BF10 ≈ 0.065，BF01 ≈ 15.5 —— 也就是说，数据其实更支持"A、B 没有差异"这个原假设，支持程度大约是 H1 的 15.5 倍。这个结果和很多人的直觉（"B 组转化率数字更高，肯定是 B 更好"）是相反的，下一节会解释为什么。

 5. 后验分析与蒙特卡洛采样

post_a_A, post_b_A = x_A + a_prior, (n_A - x_A) + b_prior
post_a_B, post_b_B = x_B + a_prior, (n_B - x_B) + b_prior

N_SAMPLES = 200_000
samples_pA = np.random.beta(post_a_A, post_b_A, size=N_SAMPLES)
samples_pB = np.random.beta(post_a_B, post_b_B, size=N_SAMPLES)
samples_delta = samples_pB - samples_pA

p_superior = np.mean(samples_delta > 0.0)   # B 优于 A 的概率
rope = 0.002
p_rope = np.mean(np.abs(samples_delta) <= rope)

结果：P(pB > pA) ≈ 90%，P(落在 ROPE 内) ≈ 15.7%，delta 的 95% 可信区间约为 [-0.31%, +1.49%]。

⚠️ 五、关键陷阱BF打架

看到这里你可能已经发现问题了：贝叶斯因子说"支持无差异"，但后验概率却说"B 有 90% 概率更好"，这两个结论完全不是一回事，必须放在一起解读，不能只挑一个说服自己。原因在于它们回答的根本是不同的问题：

• BF10 回答的是一个很苛刻的问题："pA 和 pB 完全相等"，相对于"pA、pB 各自独立、可以是任意值"，哪个更被数据支持。当真实效应很小、又要在整个概率空间上摊先验时，"完全相等"这个简单假设反而不容易被推翻——这就是统计学里经典的 Lindley/Jeffreys 悖论：样本量越大，点假设检验对"小效应"反而越保守。

• P(pB > pA) 回答的是方向问题："如果真有差异，更可能偏向哪边"，哪怕真实差异只有 0.01%，只要数据稍微偏一点，这个概率也可能冲到 80%~90%，它完全不衡量"差异到底有没有意义"。

用双比例 z 检验交叉验证一下（这里的 p-value ≈ 0.20），结论和 BF 是一致的：在传统假设检验框架下，这个差异也没有达到统计显著。

怎么用才对：不要只看某一个指标下结论，三件事一起看——BF（差异是否值得当回事）、方向概率（如果有差异，偏向哪边）、ROPE（差异是否有业务意义）。本文进一步做了一个稳健性检查：固定"真实"转化率不变（A=5%，B=6%），重复模拟 200 次数据，结果发现即便真实差异确实存在，也只有约 26% 的重复实验会让 BF10 落在支持 H1 的一侧。这说明在当前样本量和先验设定下，点假设形式的 BF 检验对这个效应量的"功效"本身就偏低——不是代码算错了，而是这本来就是这个检验方法在小效应量场景下的正常表现，值得在下结论前心里有数。

📊 六、五大可视化分析图解

完整代码如下

import numpy as np
import pandas as pd
from scipy.stats import beta as beta_dist
from scipy.special import betaln
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import seaborn as sns

# ================= 1. 中文字体与图表风格 =================
# 设置全局字体，解决中文显示问题
sns.set_theme(style="whitegrid")
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# ================= 2. 生成虚拟 A/B 数据 =================
np.random.seed(42)
p_true_A = 0.050
p_true_B = 0.060
n_A = 5000
n_B = 4800

A = np.random.binomial(1, p_true_A, size=n_A)
B = np.random.binomial(1, p_true_B, size=n_B)

x_A = A.sum()  # A 组成功数
x_B = B.sum()  # B 组成功数
print(f"观测: A 组 {x_A}/{n_A} ({x_A/n_A:.2%}), B 组 {x_B}/{n_B} ({x_B/n_B:.2%})")


# ================= 3. 贝叶斯因子函数 =================
def bayes_factor_beta_binomial(xA, nA, xB, nB, a=1, b=1):
    """默认 Beta(1,1) 无信息先验；正式分析中请显式传入业务先验，
    避免"函数默认值"和"正文实际使用值"不一致造成误解。"""
    num = betaln(xA + a, nA - xA + b) + betaln(xB + a, nB - xB + b)
    den = betaln(a, b) + betaln(xA + xB + a, (nA + nB) - (xA + xB) + b)
    return np.exp(num - den)


def beta_params_from_mean_strength(mean, strength):
    """按'先验均值 mean + 先验强度 strength(=a+b)'两个可解释的业务维度，
    反推 Beta 分布的 (a, b) 参数。strength 越大先验越"固执"。"""
    a = mean * strength
    b = (1 - mean) * strength
    return a, b


# ---- 先验设定：围绕历史基线转化率展开，而不是围绕0.5 ----
# 假设业务历史基线转化率约为 5%（可替换为真实历史数据），
# 先验强度取 4，相当于"大约4次历史观测"的信息量，是比较温和、克制的先验。
PRIOR_MEAN = 0.05
PRIOR_STRENGTH = 4
a_prior, b_prior = beta_params_from_mean_strength(PRIOR_MEAN, PRIOR_STRENGTH)
print(f"先验: Beta({a_prior:.2f}, {b_prior:.2f})，均值={PRIOR_MEAN:.1%}，强度={PRIOR_STRENGTH}")

BF10 = bayes_factor_beta_binomial(x_A, n_A, x_B, n_B, a=a_prior, b=b_prior)
print(f"贝叶斯因子 BF10 = {BF10:.4f} (log10={np.log10(BF10):.3f})")
print(f"  → BF01 = {1/BF10:.2f}，即数据对 H0(pA=pB) 的支持是 H1 的 {1/BF10:.1f} 倍")

# ================= 4. 后验分布与蒙特卡洛采样 =================
post_a_A = x_A + a_prior
post_b_A = (n_A - x_A) + b_prior
post_a_B = x_B + a_prior
post_b_B = (n_B - x_B) + b_prior

N_SAMPLES = 200_000
samples_pA = np.random.beta(post_a_A, post_b_A, size=N_SAMPLES)
samples_pB = np.random.beta(post_a_B, post_b_B, size=N_SAMPLES)
samples_delta = samples_pB - samples_pA

p_superior = np.mean(samples_delta > 0.0)
# ROPE 的选取需要业务依据，这里假设"转化率差异小于0.2个百分点，
# 对应的流量/收入变化在当前体量下不足以覆盖切换成本"，因此定义为无感区间。
# 实际使用时请替换为真实的业务测算。
ROPE = 0.002
p_rope = np.mean(np.abs(samples_delta) <= ROPE)
ci_delta = np.percentile(samples_delta, [2.5, 97.5])

print(f"P(p_B > p_A | data) = {p_superior:.3%}")
print(f"后验等价(ROPE±{ROPE:.3%})概率 = {p_rope:.3%}")
print(f"差值 delta 的 95% 可信区间 = [{ci_delta[0]:.4%}, {ci_delta[1]:.4%}]")

# ================= 5. 稳健性检查：重复模拟多次，看 BF10 有多"抖" =================
# 单次模拟的结果会受随机种子影响，尤其是在效应量小、样本量有限时。
# 这里固定"真实"转化率不变，重复抽样很多次，观察 BF10 的分布，
# 用来说明"这次抽到的样本恰好让 BF10 支持 H0"这件事本身有多大概率发生。
N_REPS = 200
bf_reps = []
for seed in range(N_REPS):
    rng = np.random.default_rng(seed)
    xa = rng.binomial(1, p_true_A, n_A).sum()
    xb = rng.binomial(1, p_true_B, n_B).sum()
    bf_reps.append(bayes_factor_beta_binomial(xa, n_A, xb, n_B, a=a_prior, b=b_prior))
bf_reps = np.array(bf_reps)
frac_favor_h1 = np.mean(bf_reps > 1)
print(f"[稳健性] 重复模拟{N_REPS}次: BF10 中位数={np.median(bf_reps):.3f}, "
      f"5%-95%分位=[{np.percentile(bf_reps,5):.3f}, {np.percentile(bf_reps,95):.3f}]")
print(f"[稳健性] 在真实差异确实存在(5% vs 6%)的前提下，"
      f"仅有 {frac_favor_h1:.1%} 的重复实验会让 BF10>1(支持H1)。")
print("  → 这说明当前样本量下，点假设(pA=pB)形式的BF检验本身对这个效应量偏保守，"
      "该结论应与后验方向概率、ROPE 一起解读，不能只看BF一个指标。")

# ================= 6. 五大可视化分析 =================

# --- 图1：先验与后验分布对比 ---
plt.figure(figsize=(10, 6))
x = np.linspace(0, 0.15, 1000)
plt.plot(x, beta_dist.pdf(x, a_prior, b_prior), 'k--',
         label=f'先验分布 Beta({a_prior:.1f},{b_prior:.1f})\n(均值={PRIOR_MEAN:.0%}, 强度={PRIOR_STRENGTH})',
         linewidth=2)
plt.plot(x, beta_dist.pdf(x, post_a_A, post_b_A), 'b-', label=f'后验分布 A组 (n={n_A})', linewidth=2)
plt.plot(x, beta_dist.pdf(x, post_a_B, post_b_B), 'r-', label=f'后验分布 B组 (n={n_B})', linewidth=2)
plt.title('图1: 先验 vs 后验分布对比', fontsize=15, fontweight='bold')
plt.xlabel('转化率 (p)', fontsize=12)
plt.ylabel('概率密度', fontsize=12)
plt.xlim(0, 0.12)
plt.legend(fontsize=11)
plt.tight_layout()
plt.savefig('fig1_prior_posterior.png', dpi=150)
plt.close()

# --- 图2：后验差值分布 ---
plt.figure(figsize=(10, 6))
sns.kdeplot(samples_delta, fill=True, color='purple', alpha=0.3)
plt.axvline(0, color='red', linestyle='--', label='无差异线 (0)')
plt.axvline(ci_delta[0], color='gray', linestyle=':', label='95%可信区间')
plt.axvline(ci_delta[1], color='gray', linestyle=':')
plt.title(f'图2: 差值分布 ($p_B - p_A$)\n'
          f'P(B优于A)={p_superior:.1%}  |  BF10={BF10:.3f}(支持H0)  |  P(ROPE内)={p_rope:.1%}',
          fontsize=14, fontweight='bold')
plt.xlabel('转化率差值', fontsize=12)
plt.ylabel('概率密度', fontsize=12)
plt.legend(fontsize=11)
plt.tight_layout()
plt.savefig('fig2_delta_distribution.png', dpi=150)
plt.close()

# --- 图3：连续监测轨迹（带证据强度背景色块） ---
plt.figure(figsize=(10, 6))
step = 100
n_steps = min(n_A, n_B) // step
bf_trajectory, sample_sizes = [], []
for i in range(1, n_steps + 1):
    curr_n = i * step
    curr_xA = A[:curr_n].sum()
    curr_xB = B[:curr_n].sum()
    bf = bayes_factor_beta_binomial(curr_xA, curr_n, curr_xB, curr_n, a_prior, b_prior)
    bf_trajectory.append(bf)
    sample_sizes.append(curr_n)

log_bf = np.log10(bf_trajectory)
y_min, y_max = min(log_bf.min(), -2.5), max(log_bf.max(), np.log10(10) + 0.3)

# 证据强度背景色块：<0 支持H0；0~log10(3) 弱证据；log10(3)~log10(10) 中等；>log10(10) 强证据
plt.axhspan(y_min, 0, color='#a6c8e0', alpha=0.25, label='支持 H0(无差异)')
plt.axhspan(0, np.log10(3), color='#f0e68c', alpha=0.3, label='证据较弱')
plt.axhspan(np.log10(3), np.log10(10), color='#f4b183', alpha=0.3, label='中等支持H1')
plt.axhspan(np.log10(10), y_max, color='#e06666', alpha=0.3, label='强支持H1')

plt.plot(sample_sizes, log_bf, 'b-', linewidth=2, zorder=5)
plt.axhline(0, color='gray', linestyle='--', linewidth=1)
plt.axhline(np.log10(3), color='orange', linestyle='--', linewidth=1)
plt.axhline(np.log10(10), color='red', linestyle='--', linewidth=1)
plt.ylim(y_min, y_max)
plt.title('图3: 连续监测轨迹（背景色块=证据强度分区）', fontsize=15, fontweight='bold')
plt.xlabel('样本量', fontsize=12)
plt.ylabel('Log10(贝叶斯因子 BF10)', fontsize=12)
plt.legend(fontsize=10, loc='upper right')
plt.tight_layout()
plt.savefig('fig3_sequential_monitoring.png', dpi=150)
plt.close()

# --- 图4：后验预测分布 ---
plt.figure(figsize=(10, 6))
future_n = 1000
pred_A = np.random.binomial(future_n, samples_pA)
pred_B = np.random.binomial(future_n, samples_pB)
pred_delta = pred_B - pred_A

sns.histplot(pred_delta, bins=50, kde=True, color='teal')
plt.title(f'图4: 后验预测 (未来 {future_n} 次展示的预期差异)', fontsize=15, fontweight='bold')
plt.xlabel('B组比A组多出的转化数', fontsize=12)
plt.ylabel('频次', fontsize=12)
plt.tight_layout()
plt.savefig('fig4_posterior_predictive.png', dpi=150)
plt.close()

# --- 图5：先验敏感性分析（固定先验均值，只变强度；对数坐标） ---
plt.figure(figsize=(10, 6))
prior_strengths = [1, 2, 4, 10, 20, 50, 100]
bf_values = []
for s in prior_strengths:
    a_temp, b_temp = beta_params_from_mean_strength(PRIOR_MEAN, s)
    bf_temp = bayes_factor_beta_binomial(x_A, n_A, x_B, n_B, a=a_temp, b=b_temp)
    bf_values.append(bf_temp)

plt.plot(prior_strengths, bf_values, 'ro-', linewidth=2, markersize=8)
plt.axhline(1, color='gray', linestyle='--', label='BF=1 (无差异)')
plt.yscale('log')
plt.title(f'图5: 先验敏感性分析（先验均值固定={PRIOR_MEAN:.0%}，仅改变先验强度）',
          fontsize=14, fontweight='bold')
plt.xlabel('先验强度 (a + b)', fontsize=12)
plt.ylabel('贝叶斯因子 (BF10, 对数坐标)', fontsize=12)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('fig5_prior_sensitivity.png', dpi=150)
plt.close()

# ================= 7. 结果汇总表 =================
summary = pd.DataFrame({
    "指标": ["A组观测转化率", "B组观测转化率", "BF10", "BF01(=1/BF10)",
             "P(pB>pA|data)", f"P(|delta|<={ROPE:.1%}) [ROPE]",
             "delta 95%可信区间", f"重复{N_REPS}次中BF10>1的比例"],
    "数值": [f"{x_A/n_A:.3%}", f"{x_B/n_B:.3%}", f"{BF10:.4f}", f"{1/BF10:.2f}",
             f"{p_superior:.3%}", f"{p_rope:.3%}",
             f"[{ci_delta[0]:.4%}, {ci_delta[1]:.4%}]", f"{frac_favor_h1:.1%}"]
})
print("\n===== 结果汇总 =====")
print(summary.to_string(index=False))

图1：先验与后验分布对比 展示数据如何让分布曲线"移动与收缩"，先验曲线不再对称地堆在 0.5 附近，而是贴着业务先验均值展开，更贴近转化率的真实取值范围。B 组的后验整体略右移（转化率更高），且区间变窄（不确定性降低）。

图2：后验差值分布（δ = pB − pA） 直接看差值的密度图，标题里同时给出 BF10、方向概率、ROPE 概率三个指标，避免读者只盯着"B 优于 A 的概率"一个数字下结论。

图3：连续监测轨迹（Sequential Monitoring） 模拟"边跑边看"的场景。X 轴是样本量，Y 轴是 log10(BF)。背景色块按证据强度分区（支持 H0 / 证据弱 / 中等支持 H1 / 强支持 H1），曲线落在哪个色块一目了然。贝叶斯因子相对更适合连续监测，但也并非完全没有"偷看"风险，严谨场景建议结合预注册的先验和分析计划。

图4：后验预测分布（Posterior Predictive） 预测"未来 1000 次展示会有多少转化差异"，从"现在谁好"进一步延伸到"接下来会发生什么"，辅助业务做资源预案。

图5：先验敏感性分析 与原来"围绕 0.5 改变强度"不同，这里固定先验均值在业务基线（5%），只改变先验强度，纵轴改为对数坐标，能更清楚地看到：先验越"固执"（强度越大），BF10 越向"无差异"方向收缩——这本身也是需要在报告里说明的一个稳健性维度。

📌 七、总结

贝叶斯检验不仅仅是统计学方法，更是一种贴近现实的决策思维：

1. 带着常识看证据（先验要贴合业务实际，不要图省事围绕 0.5 展开）。

2. 用贝叶斯因子量化支持力度，但要清楚它回答的是"要不要认为有差异"，不是"哪边更好"。

3. 用后验分布刻画不确定性（方向概率、ROPE、可信区间，三件套一起看）。

4. 对单次模拟结果保持警惕，必要时用重复模拟检查结论的稳健性。

5. 用后验预测直接回答未来，把统计结论落到业务决策上。