Pyro用PyTorch玩转概率编程与贝叶斯推断

可视化演示

贝叶斯思想

老主厨收徒，从不先教菜谱，而是先教“试菜哲学”：你心里对一道菜有个预期（先验），每让客人尝一口就获得一条反馈（似然），把预期和反馈揉在一起，才形成你对这道菜真正的把握（后验）。贝叶斯推断的全部奥义就在这个循环里。

在概率语言中，假设我们要估计的参数是 \(\theta\)（比如盐的克数）。我们先给它一个先验分布 \(p(\theta)\)，它代表“旧配方”的信念。看到数据 \(D\) 后，通过似然 \(p(D|\theta)\) 刻画“在某个配方下客人给出这种反馈的概率”。贝叶斯定理把二者合成后验：

厨师的故事告诉我们一个深刻结论：后验不是“一个最佳配方”，而是一整张“配方分布图”。SVI 是快速逼近这张图的学徒，MCMC 是仔细品尝每个角落的老师傅。

---

Pyro初体验

厨房里来了个新助手叫 Pyro。它不直接做菜，而是帮你搭“菜谱的骨架”。在概率编程里，你只需要用几行代码描述“数据是怎么生成的”，剩下的推断交给框架。

Pyro 的模型就是一个普通的 Python 函数。函数体内出现两个核心原语：pyro.sample 表示“从一个分布里抽一个随机变量”，pyro.observe 表示“这个随机变量有观测值，要拿它对比数据”。在这两个原语之间，PyTorch 的张量运算照常进行。

1. sample 声明潜变量：torch.Size([]) 表示标量，dist.Normal(3.5, 1.0) 是标准正态分布，也就是“旧配方”。

$$\theta \sim \text{Normal}(3.5, 1.0)$$

2. plate 批量向量化：plate 将独立同分布的数据维度声明出来，让 Pyro 能做向量化批处理，等价于假设样本间条件独立。

$$\prod_{i=1}^n p(y_i \mid \theta)$$

3. observe 绑定观测：obs=data 将随机变量 obs 固定为真实数据，使 Pyro 在推断时自动应用贝叶斯定理。

$$p(\theta \mid D) \propto p(\theta) \prod_{i=1}^n p(y_i \mid \theta)$$

Pyro 的优雅在于：你只写“生成故事”，不写“推断算法”。推断是助手的事——你指定用 SVI 还是 MCMC，它就把后验端上来。

---

def model(data):    theta = pyro.sample("theta", dist.Normal(3.5, 1.0))    with pyro.plate("data", len(data)):        pyro.sample("obs", dist.Normal(theta, 1.0), obs=data)

PyTorch地基

Pyro 能如此灵活，功劳首先归功于 PyTorch 这个“万用食材柜”。张量就是你的食材盘——无论是一维向量还是高维矩阵，都能整齐排列、高速运算。自动微分则是“自动调味机”：你只需写清楚前向过程，反向传播自动帮你算出每个参数的偏导。

Pyro 在 PyTorch 之上加了一层“概率调料架”：torch.distributions 提供了几十种分布（正态、泊松、Dirichlet 等），这些分布既可以作为 Pyro 的 dist，也能单独用于张量采样。

| PyTorch 能力 | 在 Pyro 中的角色 |

|:——|:——|

| 张量（Tensor） | 存储参数、数据和采样结果 |

| 自动微分（autograd） | SVI 的梯度下降自动求导 |

| 分布模块（distributions） | 构建先验、似然与 guide 分布 |

| DataLoader | 大数据下分批次喂给模型 |

1. 张量统一数据形态：所有参数和观测都化为 torch.Tensor，支持 GPU 并行。

$$\mathbf{X} \in \mathbb{R}^{n \times d}$$

2. 自动微分计算梯度：ELBO 对变分参数 \(\phi\) 求导，反向传播自动完成。

$$\nabla_{\phi} \mathcal{L}(\phi)$$

3. 分布模块统一接口：每个分布都有 sample 和 log_prob，便于 Pyro 统一处理随机变量。

$$q_{\phi}(\theta) \to \text{sample}, \log q_{\phi}(\theta)$$

这一层地基决定了 Pyro 的扩展性：当你把贝叶斯层叠加到深度网络上时，PyTorch 的模块化让“深度概率模型”就像搭积木。

---

贝叶斯回归

回到厨房：你想搞明白“火候”（\(x\)）和“客人评分”（\(y\)）之间的关系。经典回归只给一条最佳拟合线，而贝叶斯回归给出“线的分布”——每条可能的线都有它的可信度。

贝叶斯线性回归假设：

对权重 \(w\) 和偏置 \(b\) 设置先验 \(w \sim \text{Normal}(0, 1)\)，然后让数据说话。

相比最小二乘，贝叶斯回归的产出是后验分布。你不仅能说“斜率约为 3.2”，还能说“斜率落在 [2.8, 3.6] 的可信度是 95%”——这在风控和医疗场景里意味着决策的安全性。

---

SVI变分推断

贝叶斯推断的难点在于分母积分。当参数维度高时，精确计算是不可能的。SVI（随机变分推断）的解决思路是“雇佣一位快速学徒”：用一族简单的分布 \(q_{\phi}(\theta)\)（guide）去近似真实后验 \(p(\theta|D)\)。学徒在一次次对照中调整自己的动作参数 \(\phi\)，让 \(q\) 越来越像 \(p\)。

衡量两个分布差距的工具是 KL 散度。直接最小化 KL 不可行，所以我们最大化它的等价下界 ELBO：

SVI 的优势是快——哪怕模型有百万参数，它也能在 GPU 上用梯度下降逼近后验。代价是近似，因为 \(q\) 的分布族可能不够灵活。

---

MCMC采样

如果说 SVI 是快速学徒，MCMC 就是老师傅。它不假设分布形状，而是直接在后验分布上“游走采样”：每走一步就留下一个样本，走过足够多步之后，留下的样本集合就近似后验分布本身。

Pyro 默认推荐 NUTS 内核。它的巧妙之处在于自动决定“每一步跳多远”：在参数空间中按照哈密顿动力学优雅地滑动，避免普通随机游走那种原地打转的低效。

1. 马尔可夫链：下一个样本只依赖当前样本，构成一条链。

$$\theta^{(t+1)} \sim T(\theta^{(t)} \to \theta^{(t+1)})$$

2. 平稳分布：设计转移核 \(T\) 使链的平稳分布恰好是后验。

$$\pi(\theta) = p(\theta \mid D)$$

3. NUTS 自动步长：利用梯度信息在相空间模拟动力学，避免随机游走低效问题。

$$\theta \leftarrow \text{Leapfrog}(\theta, r, \epsilon)$$

4. 采样平均：舍弃前 \(B\) 个老化样本后，用经验分布近似后验。

$$\hat{p}(\theta \mid D) \approx \frac{1}{T}\sum_{t=B+1}^{T}\delta(\theta - \theta^{(t)})$$

MCMC 的代价是慢，但换来的是渐近精确。当后验分布多峰、扭曲、非标准时，老师傅依然能品尝出每个角落的味道。

---

预测与前沿

学好了菜谱，最终要用它开新店。在贝叶斯世界里，对新输入的预测不是点估计，而是一个后验预测分布：

在 Pyro 中，预测和前沿应用都站在同一块地基上：后验是食材，分布是刀具，而 do 算子是分子料理的喷雾器——它帮你回答“如果我把盐换糖，客人评分会怎样”这类反事实问题。

---

AI 应用场景

Python 演示

import torchimport pyroimport pyro.distributions as distfrom pyro.infer import SVI, Trace_ELBO, MCMC, NUTS, Predictive pyro.set_rng_seed(0) # 生成模拟数据：真实关系 y = 3.2x + 1.8 + 噪声x = torch.linspace(0, 10, 50)y = 3.2 * x + 1.8 + torch.randn(50) * 1.5  def model(x, y):    w = pyro.sample("w", dist.Normal(0.0, 5.0))    b = pyro.sample("b", dist.Normal(0.0, 5.0))    sigma = pyro.sample("sigma", dist.HalfCauchy(1.0))    mu = w * x + b    with pyro.plate("data", len(x)):        pyro.sample("obs", dist.Normal(mu, sigma), obs=y)  def guide(x, y):    w_mu = pyro.param("w_mu", torch.tensor(1.0))    w_sigma = pyro.param("w_sigma", torch.tensor(1.0),                         constraint=dist.constraints.positive)    b_mu = pyro.param("b_mu", torch.tensor(1.0))    b_sigma = pyro.param("b_sigma", torch.tensor(1.0),                         constraint=dist.constraints.positive)    sigma = pyro.param("sigma_q", torch.tensor(1.0),                       constraint=dist.constraints.positive)    pyro.sample("w", dist.Normal(w_mu, w_sigma))    pyro.sample("b", dist.Normal(b_mu, b_sigma))    pyro.sample("sigma", dist.LogNormal(torch.tensor(0.0), sigma))  def main():    # SVI 变分推断    svi = SVI(model, guide, pyro.optim.Adam({"lr": 0.05}), loss=Trace_ELBO())    for step in range(1000):        svi.step(x, y)    print("SVI w_mu =", pyro.get_param_store()["w_mu"].item())     # MCMC 采样(NUTS内核)    nuts = NUTS(model)    mcmc = MCMC(nuts, num_samples=200, warmup_steps=100)    mcmc.run(x, y)    samples = mcmc.get_samples()    print("MCMC w mean =", samples["w"].mean().item())     # 后验预测    predictive = Predictive(model, samples, return_sites=("obs",))    x_new = torch.tensor([11.0, 12.0])    pred = predictive(x_new, None)    print("预测均值:", pred["obs"].mean(dim=0).tolist())    print("预测标准差:", pred["obs"].std(dim=0).tolist())  if __name__ == "__main__":    main()

1. pyro.sample("w", dist.Normal(0.0, 5.0))将权重 w 声明为随机变量，先验设为均值为 0 的正态分布。
2. with pyro.plate("data", len(x))批量处理所有数据点，声明它们条件独立，提升向量化效率。
3. pyro.param(...constraint=positive)声明变分参数时用正值约束，保证标准差合法。
4. SVI(model, guide, optimizer, loss=Trace_ELBO())组装学徒（guide）、考试题（ELBO）和优化器，开始训练。
5. MCMC(NUTS(model), num_samples=200, warmup_steps=100)请出老师傅，先烧 100 步老火，再取 200 个后验样本。
6. Predictive(model, samples, return_sites=("obs",))用后验样本做预测，输出的 obs 就是后验预测分布。

常见误区

面试题