贝叶斯推断与概率图模型:从共轭先验到 MCMC 实战

贝叶斯推断把参数视为随机变量,由先验与似然共同决定后验分布。本文讲清共轭先验(Beta-Binomial、Normal-Normal、Dirichlet)、MLE 与 MAP 的关系、后验预测分布、MCMC 采样(Metropolis-Hastings、Gibbs、NUTS)与 PyMC 实战,并覆盖贝叶斯网络、条件独立与 d-分离、朴素贝叶斯、变分推断 ELBO 及贝叶斯 A/B 测试落地。

引言

频率学派把参数看作固定的未知常数,用重复采样的频率解释概率;贝叶斯学派则把参数本身看作随机变量,用先验分布表达在看到数据前的信念,用似然吸收数据证据,得到后验分布。后验不只是给出一个点估计,而是给出参数的完整不确定性刻画——这在样本少、需要决策风险量化的场景里价值巨大。

代价是后验往往没有解析解,必须借助采样或近似推断。本文从贝叶斯定理出发,讲清共轭先验的便利、MCMC 与变分推断两大推断路线,以及概率图模型如何用图结构编码条件独立性,最后用 PyMC 完成端到端实战。

前置:概率分布与假设检验见 /ml-model-evaluation/;因果图与 do 算子的区别见 /ml-causal-inference/。

目录

1. 贝叶斯定理与三要素

贝叶斯定理把联合概率的两种分解联系起来:

P(θ | D) = P(D | θ) · P(θ) / P(D)

四个部分各有名字:

名称符号含义
先验P(θ)看到数据前对参数的信念
似然P(D | θ)给定参数时观测到数据的概率
后验P(θ | D)看到数据后更新的信念
证据P(D)归一化常数,∫P(D|θ)P(θ)dθ

分母 P(D) 与 θ 无关,常写作正比形式:

后验  ∝  似然 × 先验
import numpy as np
from scipy import stats

# 抛硬币:先验 Beta(1,1)(均匀),观测 7 次正面 3 次反面
prior_a, prior_b = 1, 1
heads, tails = 7, 3
post_a, post_b = prior_a + heads, prior_b + tails   # Beta 共轭更新
x = np.linspace(0, 1, 100)
print("后验均值:", round(post_a / (post_a + post_b), 4))
print("后验 95% 区间:", stats.beta.ppf([0.025, 0.975], post_a, post_b).round(3))

2. 共轭先验

当先验与后验属于同一分布族时,称先验为共轭先验(Conjugate Prior),此时后验有闭式解,无需采样。这是贝叶斯推断里最省事的一类情形。

似然共轭先验后验更新
Bernoulli / BinomialBeta(α, β)α+成功次数, β+失败次数
PoissonGamma(α, β)α+Σx, β+n
Normal(方差已知)Normal(μ₀, σ₀²)精度加权平均
Normal(均值已知)Inverse-Gamma形状、尺度累加
Categorical / MultinomialDirichlet(α)α_k + 类别计数

2.1 Beta-Binomial 的直觉

Beta(α, β) 可以理解成「先验里已经看过 α-1 次成功、β-1 次失败」。观测数据直接加到参数上,先验越弱(α, β 越小),数据影响力越大。

def beta_update(alpha, beta, successes, failures):
    return alpha + successes, beta + failures

# 强先验 Beta(50,50),观测 7/3
print(beta_update(50, 50, 7, 3))    # (57, 53) —— 数据被先验拉回
# 弱先验 Beta(1,1)
print(beta_update(1, 1, 7, 3))      # (8, 4)  —— 数据主导

2.2 Dirichlet-Multinomial

多类别场景下 Dirichlet 是 Beta 的推广:

import numpy as np

alpha = np.array([1.0, 1.0, 1.0])       # 三类先验
counts = np.array([10, 5, 2])            # 观测计数
posterior = alpha + counts
print("后验参数:", posterior)
print("后验均值:", (posterior / posterior.sum()).round(3))

3. MLE、MAP 与后验预测

3.1 三种点估计的关系

方法目标与先验关系
MLEmax P(D | θ)不使用先验
MAPmax P(θ | D) ∝ P(D|θ)P(θ)用先验作为正则
后验均值E[θ | D]用整个后验

MAP 估计等价于「似然 + 先验正则项」。例如正态先验下 MAP 就是 L2 正则(岭回归),拉普拉斯先验下就是 L1 正则(Lasso)——正则化本质上是引入了先验信念。

from scipy.optimize import minimize

def neg_log_posterior(theta, data, prior_mu=0.0, prior_sigma=1.0):
    n = len(data)
    ll = -0.5 * np.sum((data - theta) ** 2)              # 高斯似然
    lp = -0.5 * ((theta - prior_mu) / prior_sigma) ** 2  # 高斯先验
    return -(ll + lp)

data = np.array([1.2, 0.8, 1.5, 0.9, 1.1])
res = minimize(neg_log_posterior, x0=[0.0], args=(data,))
print("MAP 估计:", round(res.x[0], 4))
print("MLE(样本均值):", round(data.mean(), 4))

先验均值 0 把 MAP 估计往 0 拉,因此 MAP 小于 MLE——这就是收缩(Shrinkage)。

3.2 后验预测分布

贝叶斯推断的最终目标是预测,而不是估计参数。后验预测分布把参数不确定性积分掉:

P(x_new | D) = ∫ P(x_new | θ) P(θ | D) dθ

比「用点估计代入」更稳健,因为它对参数的不确定性做了平均,避免过拟合。

4. 概率图模型

概率图模型(Probabilistic Graphical Model, PGM)用图结构紧凑地表示多变量联合分布,节点是随机变量,边是依赖关系。

4.1 有向图:贝叶斯网络

贝叶斯网络用有向无环图(DAG)表示因果/依赖方向,联合分布按父节点分解:

P(X1, ..., Xn) = Π P(Xi | Parents(Xi))

例如「下雨 → 洒水器 → 草地湿」的网络:

from pgmpy.models import DiscreteBayesianNetwork
from pgmpy.factors.discrete import TabularCPD

model = DiscreteBayesianNetwork([("Rain", "Wet"), ("Sprinkler", "Wet")])
cpd_rain = TabularCPD("Rain", 2, [[0.8], [0.2]])
cpd_spr = TabularCPD("Sprinkler", 2, [[0.5], [0.5]])
cpd_wet = TabularCPD("Wet", 2,
                     [[1.0, 0.1, 0.1, 0.01], [0.0, 0.9, 0.9, 0.99]],
                     evidence=["Rain", "Sprinkler"], evidence_card=[2, 2])
model.add_cpds(cpd_rain, cpd_spr, cpd_wet)
print("模型合法:", model.check_model())

4.2 条件独立与 d-分离

图结构编码了条件独立关系。d-分离(d-separation) 给出判定规则:若给定观测集 Z 后 X 与 Y 之间所有路径都被阻断,则 X ⊥ Y | Z。三种基本结构:

链式    X → M → Y      给定 M 则阻断
分叉    X ← M → Y      给定 M 则阻断
对撞    X → M ← Y      给定 M 反而打通(对撞偏差)

对撞结构是贝叶斯网络最容易出错的地方:两个本独立的变量,在给定共同后代后变得相关。

4.3 无向图与朴素贝叶斯

马尔可夫随机场(MRF)用无向图表示对称依赖,用势函数定义联合分布。朴素贝叶斯则是最简单的贝叶斯网络:假设特征在类别给定下条件独立:

P(y | x1..xn) ∝ P(y) Π P(xi | y)
from sklearn.naive_bayes import GaussianNB
from sklearn.datasets import load_iris
from sklearn.model_selection import cross_val_score

X, y = load_iris(return_X_y=True)
print("朴素贝叶斯 CV:", cross_val_score(GaussianNB(), X, y, cv=5).mean().round(4))

5. 蒙特卡洛与 MCMC

绝大多数真实模型的后验无解析解,需要采样。蒙特卡洛的核心思想是:只要能从未知分布采样,就能用样本均值近似期望。

5.1 蒙特卡洛积分

# 用采样估计 E[X] 与 P(X > 2),X ~ 标准正态
samples = np.random.standard_normal(100000)
print("均值估计:", round(samples.mean(), 4))
print("P(X>2) 估计:", round((samples > 2).mean(), 4))
print("P(X>2) 理论:", round(1 - stats.norm.cdf(2), 4))

5.2 Metropolis-Hastings

当无法直接采样时,构造马尔可夫链,使其平稳分布恰好是目标后验。MH 算法:

1. 从提议分布 q(θ' | θ) 采样候选 θ'
2. 计算接受率 α = min(1, [p(θ')q(θ|θ')] / [p(θ)q(θ'|θ)])
3. 以概率 α 接受,否则留在原地
def metropolis_hastings(log_pdf, init, n_samples=20000, step=0.5):
    samples = np.zeros(n_samples)
    current = init
    for i in range(n_samples):
        proposal = current + np.random.normal(0, step)
        log_ratio = log_pdf(proposal) - log_pdf(current)
        if np.log(np.random.uniform()) < log_ratio:
            current = proposal
        samples[i] = current
    return samples

# 目标:标准正态(未归一化)
log_target = lambda x: -0.5 * x ** 2
chain = metropolis_hastings(log_target, 0.0)
print("后验均值:", round(chain[5000:].mean(), 4))   # 丢弃预热

5.3 Gibbs 与 HMC/NUTS

  • Gibbs 采样:逐个变量、用其全条件分布采样,适合共轭结构多的高维模型。
  • HMC(Hamiltonian Monte Carlo):借梯度信息在参数空间「弹跳」,在高维下比随机游走的 MH 高效得多。
  • NUTS(No-U-Turn Sampler):HMC 的自适应变体,自动调步长与轨迹长度,是 PyMC/Stan 的默认采样器。

5.4 收敛诊断

MCMC 必须检查收敛,不能直接信样本:

指标含义合格阈值
R-hat链间与链内方差比< 1.01
ESS有效样本量越大越好,> 400
Trace 图链的轨迹应像毛毛虫,不漂移
自相关样本相关性迅速衰减

6. 变分推断

采样在高维大数据下仍慢。变分推断(Variational Inference, VI) 用一个简单分布 q(θ) 近似后验,把推断变成优化:

目标:最小化 KL(q(θ) ‖ p(θ|D))
等价于:最大化 ELBO = E_q[log p(D, θ)] - E_q[log q(θ)]

ELBO 分解为「期望似然」减去「KL 散度」:前者鼓励拟合数据,后者鼓励 q 接近先验。

import numpy as np

def elbo(mu_q, log_sigma_q, data, prior_sigma=1.0):
    sigma_q = np.exp(log_sigma_q)
    # 期望对数似然(高斯模型)
    exp_ll = -0.5 * np.sum((data - mu_q) ** 2) - len(data) * log_sigma_q
    # KL(q || prior)
    kl = 0.5 * (sigma_q**2 / prior_sigma**2 + mu_q**2 / prior_sigma**2
                - 1 - 2 * log_sigma_q + 2 * np.log(prior_sigma))
    return exp_ll - kl

data = np.array([1.2, 0.8, 1.5, 0.9, 1.1])
print("ELBO@mu=1.1:", round(elbo(1.1, np.log(0.4), data), 3))

VI 比 MCMC 快几个数量级,适合大规模,代价是近似可能低估后验方差。常用方法有均场近似、随机变分推断、黑盒变分推断(BBVI)。

7. PyMC 实战

用 PyMC 做贝叶斯线性回归,把「广告投入预测销量」写成概率模型。

import numpy as np
import pymc as pm

rng = np.random.default_rng(0)
n = 100
x = rng.normal(0, 1, n)
y = 2.5 * x + 1.0 + rng.normal(0, 0.5, n)

with pm.Model() as model:
    intercept = pm.Normal("intercept", mu=0, sigma=5)
    slope = pm.Normal("slope", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=1)
    mu = intercept + slope * x
    pm.Normal("y", mu=mu, sigma=sigma, observed=y)
    trace = pm.sample(draws=1000, tune=1000, chains=4, cores=4,
                      random_seed=1, progressbar=False)

print(pm.summary(trace, var_names=["intercept", "slope", "sigma"]))

输出会给出每个参数的后验均值、标准差、HDI(最高密度区间)与 R-hat。相比最小二乘的点估计,贝叶斯回归额外告诉你斜率有多确定——HDI 是否跨过 0 直接回答了「这个效应是否可信」。

7.1 后验预测检查

with model:
    post_pred = pm.sample_posterior_predictive(trace, random_seed=1,
                                               progressbar=False)
pred_mean = post_pred.posterior_predictive["y"].mean(dim=["chain", "draw"])
print("预测 RMSE:", round(np.sqrt(((pred_mean - y) ** 2).mean()), 4))

8. 贝叶斯 A/B 测试与决策

贝叶斯方法在 A/B 测试中的优势:直接给出「B 优于 A 的概率」,无需 p 值与多重比较校正,且支持序贯查看。

from scipy import stats
import numpy as np

# A: 100/1000, B: 130/1000
a_post = stats.beta(1 + 100, 1 + 900)
b_post = stats.beta(1 + 130, 1 + 870)
samples = np.random.default_rng(0).beta
sa = a_post.rvs(200000, random_state=0)
sb = b_post.rvs(200000, random_state=0)
print("P(B > A):", round((sb > sa).mean(), 4))
print("预期提升:", round((sb - sa).mean() / sa.mean() * 100, 2), "%")
print("提升 95% 区间:", np.percentile(sb - sa, [2.5, 97.5]).round(4))

8.1 损失函数决策

光看「B 更好」不够,还要看选错的期望损失:

loss_choose_b = np.maximum(sa - sb, 0).mean()
loss_choose_a = np.maximum(sb - sa, 0).mean()
print("选 B 的期望损失:", round(loss_choose_b, 5))
print("选 A 的期望损失:", round(loss_choose_a, 5))

当「选 A 的期望损失」低于业务可接受阈值时即可停止实验上线 B。

8.2 与贝叶斯优化的联系

贝叶斯优化用高斯过程代理目标函数,再用采集函数(EI/UCB)选下一个评估点——它的「后验更新」正是本文的贝叶斯推断框架,具体实现见 /ml-automl-hpo/。

9. 常见坑与选型

现象根因处理
R-hat 远大于 1.01链未收敛增加 tune、重参数化、检查模型
采样极慢高维、共轭少换 NUTS、中心化参数、用 VI
后验被先验主导先验太强或数据太少用弱信息先验,或收集更多数据
后验区间过窄VI 低估方差用全秩 VI 或换 MCMC 验证
结果不可复现随机种子未固定固定 random_seed,报告多链结果

9.1 采样还是变分

数据量小、维度低、要精确后验?
  ├── 是 → MCMC(NUTS)
  └── 否
        ├── 大数据、要快速迭代 → 变分推断 / 随机 VI
        ├── 有共轭结构 → Gibbs / 闭式解
        └── 只是要个点估计 + 正则 → MAP(等价于带正则的 MLE)

10. 总结

10.1 核心要点

  • 贝叶斯推断 = 先验 + 似然 → 后验,后验给出完整的不确定性刻画。
  • 共轭先验让后验有闭式解;MAP 等价于带先验正则的 MLE。
  • 概率图模型用图结构编码条件独立,对撞结构是最易踩的坑。
  • 无解析解时用 MCMC(精确但慢)或 变分推断(快但近似)。

10.2 何时选贝叶斯

当样本稀少、需要不确定性量化、需要序贯更新(如在线 A/B)、或需要把领域知识写进先验时,贝叶斯方法明显优于频率方法。反之,纯预测任务、样本充足时,频率学派方法往往更简单够用。想进一步理解采样与优化的数值方法,可延伸阅读 Python 专题 ,若关注算法与工程实现深度,可参考 AI/ML 专题 。

继续阅读

探索更多技术文章

浏览归档,发现更多关于系统设计、工具链和工程实践的内容。

全部文章 返回首页

「ml」更多文章

  1. 高斯过程回归与分类:不确定性建模实战
  2. MLOps 全生命周期实践:从实验到生产闭环
  3. 公平性与偏差缓解:从度量到去偏实战