前面的模型都给你一个点估计——「答案是 0.87」。但 0.87 是「我非常确定」还是「我瞎猜的」？普通神经网络分不清。贝叶斯方法的核心野心：不输出一个数，而是输出一个分布——「大概 0.87，但也可能在 0.6~0.95 之间」。今天讲清这套思想的两大计算引擎（MCMC 采样、变分推断）、它给深度学习带来的能力（不确定性量化），以及把这套数学封装成工具的概率编程语言。
你想知道一个海量数据库里某个聚合分布的形状，但没法全表扫描（积分算不出来）。MCMC 派出一个随机游走的爬虫在状态空间里走动，规则是：概率密度高的地方多停留、低的地方少停留。走够久之后，统计爬虫的「访问频率直方图」，就近似了你要的那个分布。它和你熟悉的蒙特卡洛抽样估 p99 是同一个思路——用样本代替全量。
贝叶斯推断的核心是贝叶斯公式：后验 ∝ 似然 × 先验。写成数学是 p(θ|D) = p(D|θ)·p(θ) / p(D)。分子好算（你的模型直接给），要命的是分母 p(D)——它是对所有可能参数 θ 的积分 ∫ p(D|θ)p(θ)dθ，在高维下根本算不出来（这叫「归一化常数」难题）。
MCMC 的妙处：它不需要算那个分母——每一步只看「新旧位置的概率密度比值」，分母在比值里被约掉了。最经典的 Metropolis-Hastings 算法只有四步：
朴素随机抖动在高维里效率极低（像喝醉的人乱走）。现代采样器用哈密顿蒙特卡洛（HMC）——借物理里「小球在能量曲面上滚」的思想，让提议顺着梯度走，一步迈得远又少被拒绝。但 HMC 要手调「滚多少步」，调不好就废。NUTS（No-U-Turn Sampler，Hoffman & Gelman 2011）让算法自动决定何时该停（小球开始往回掉头时就停），这就是 PyMC / Stan 今天默认采样器的由来。
import numpy as np # 从一个未归一化的目标分布采样（双峰），演示 MCMC 不需要分母 def target(x): # 只需密度的"形状"，不需积分常数 return np.exp(-(x-2)**2/0.5) + np.exp(-(x+2)**2/0.5) samples, x = [], 0.0 for _ in range(50000): x_new = x + np.random.normal(0, 1) # ② 在附近提议 a = min(1, target(x_new) / target(x)) # ③ 比值，常数约掉 if np.random.rand() < a: # ④ 掷骰子决定是否接受 x = x_new samples.append(x) # 丢掉前段"预热(burn-in)"未收敛的样本，剩下的就近似目标分布 print(np.mean(samples[5000:]), np.std(samples[5000:]))
MCMC 慢且难扩展——百万样本对大数据集是灾难。变分推断（VI）换思路：与其精确采样那个复杂后验，不如找一个简单分布（比如高斯）去「拟合」它。这就像用一个预计算的物化视图 / 缓存去近似一个昂贵的实时聚合查询——牺牲一点精度，换巨大的速度。VI 把「推断问题」变成了你最熟的「优化问题」：调参数让近似分布尽量贴近真后验。
设真后验是 p(θ|D)（算不出），引入一个带参数 φ 的简单分布 qφ(θ)（比如均值方差待定的高斯）。目标：调 φ，让 q 和 p 的「距离」最小。距离用 KL 散度衡量（Day 33 讲过：两个分布的信息差）。
但 KL(q‖p) 里又藏着那个算不出的后验。数学上一变形，最小化 KL 等价于最大化一个叫 ELBO（证据下界）的量：
直觉拆解：第一项逼 q 解释好观测数据（似然要高）；第二项是个「拉回先验」的弹簧，防止过拟合。两项拔河——这正是「拟合 vs 正则」的贝叶斯版本，你在 Day 10 见过同样的张力。ELBO 之所以叫「下界」，是因为它永远 ≤ 真实证据 log p(D)，把它顶高，就等于把 q 推向 p。
关键是 ELBO 可以求梯度，于是能用 Adam 直接跑。把高斯的「采样」改写成「均值 + 标准差 × 标准噪声」，梯度就能穿过随机采样回流——这个重参数化技巧（reparameterization trick）是 Kingma & Welling 的 VAE 论文（2013）的核心，也是 VAE、扩散模型背后的同一块积木。PyMC / NumPyro 里这套自动化版本叫 ADVI。
import pymc as pm import numpy as np y = np.random.normal(5, 2, size=200) # 假数据：真均值 5 with pm.Model() as model: mu = pm.Normal("mu", 0, 10) # 先验：均值不知道 sigma = pm.HalfNormal("sigma", 5) pm.Normal("obs", mu, sigma, observed=y) # 似然 # 不用 MCMC 采样，改用变分推断拟合 —— 大数据集快几个量级 approx = pm.fit(20000, method="advi") # 最大化 ELBO idata = approx.sample(1000) # 从拟合好的 q 里抽样 print(idata.posterior["mu"].mean().item()) # ≈ 5，且带后验区间
普通神经网络永远「自信满满」——哪怕喂它一张从没见过的乱码图，它也会斩钉截铁说「这是猫，99%」。这就像一个从不返回错误码、永远返回 200 OK 的服务，哪怕后端早就崩了。不确定性量化（UQ）给模型装上一个诚实的「我不知道」信号——区分「我见过类似数据，有把握」和「这超出我的认知范围，别信我」。
关键是区分两种本质不同的不确定性，这个区分是整个领域的地基：
怎么让神经网络给出 epistemic 信号？核心思路：别只信一个模型，要看「一群模型」的分歧。数据密集处大家答案一致（低不确定）；没见过的区域各模型乱猜、互相打架（高不确定）。两种最实用的实现：
import torch # MC Dropout：推理时保持 dropout 打开，跑多次拿分歧 def predict_with_uncertainty(model, x, n=30): model.train() # 关键：train 模式让 dropout 仍激活 preds = torch.stack([model(x) for _ in range(n)]) # N 次随机前向 mean = preds.mean(0) # 预测值 std = preds.std(0) # 不确定性：分歧大 = 模型没把握 return mean, std mean, std = predict_with_uncertainty(my_net, x_test) # std 高的样本 → 模型在"瞎猜"，应转人工审核或拒答 mask = std.squeeze() > 0.3 print(f"需人工复核的样本数: {mask.sum().item()}")
前面三个概念的数学都很硬。概率编程语言（PPL）是把它们封装成一门声明式语言：你只「声明」概率模型（先验是什么、数据怎么生成），推断引擎自动帮你跑 MCMC 或 VI。这正是 SQL 之于数据库 的关系——你写「想要什么（what）」，查询优化器决定「怎么算（how）」。PyMC / Stan 就是概率世界的 SQL，把你从手写采样器里解放出来。
没有 PPL 之前，每换一个模型你都得重写 Metropolis 接受率、梯度、收敛诊断——极易出错。PPL 的核心机制是把模型定义和推断算法解耦：
主流选择：Stan——统计学界金标准、最严谨，独立 DSL；PyMC——纯 Python、上手最快；NumPyro——架在 JAX 上、GPU 加速、最快，适合大模型贝叶斯。它们共享同一套引擎核心：自动微分 + NUTS 采样 + 变分推断，区别只在语法与速度。
import pymc as pm import numpy as np # 贝叶斯线性回归：注意我们只"声明"模型，不写任何采样逻辑 x = np.linspace(0, 10, 50) y = 2.5 * x + 1.0 + np.random.normal(0, 1, 50) # 真斜率 2.5 with pm.Model() as model: slope = pm.Normal("slope", 0, 5) # ① 先验 inter = pm.Normal("inter", 0, 5) noise = pm.HalfNormal("noise", 2) # ② 似然 + ③ 绑数据 pm.Normal("y", slope*x + inter, noise, observed=y) idata = pm.sample(1000) # ④ 引擎自动 NUTS 采样 + 诊断 # 得到的不是一个斜率，而是斜率的整个后验分布 pm.summary(idata) # 输出均值、94% 区间、r_hat、ess