贝叶斯校准 MCMC 模拟器

Q: 为什么需要 burn-in？应丢弃多少步？

MCMC 的初值往往远离后验高密度区，最初的几百到几千步是「向高密度区漂移」的过渡样本。如果把它们混入统计，会让后验均值与方差产生偏差，因此需要丢弃，这就是 burn-in。本工具默认 500 步对于弱信息先验与数千样本的链子是较保守的取值。观察 trace plot，确认左端的漂移已消失即可。严谨场合需要并行多条链并用 Gelman-Rubin 的 R̂ < 1.01 判断收敛。

Q: 有效样本数 ESS 偏小为什么是个问题？

MCMC 样本之间存在自相关，因此 5000 个原始样本「折算为独立样本数」可能只有几百甚至几十，这就是 ESS。后验均值的标准误为 σ_post/√ESS，当 ESS 跌破 100 时，估计值的标准误就超过后验标准差的 10%，结果的小数点第二位以下就不能信任了。提高 ESS 的三种途径是：(1) 优化提案宽度；(2) 增加样本数；(3) 切换到 Hamiltonian Monte Carlo / NUTS 等更高效的采样器。

参数设置

先验均值 μ_p

参数 θ 的先验知识中心

先验标准差 σ_p

越大越弱信息（由数据主导）

观测噪声 σ_y

观测模型的已知标准差

观测均值 ȳ

N 个观测的样本均值

观测数 N

后验 std 按 1/√N 收缩

MH 提案宽度 step

约 2.38·σ_post 时最优

MCMC 样本数

Burn-in

初始过渡样本的丢弃数

计算结果

—

后验均值 μ_post

—

后验标准差 σ_post

—

95% 可信区间宽度

—

MH 接受率 (%)

—

有效样本数 ESS

—

信息增益 (nats)

—

先验·似然·后验叠加 + MH walker

蓝色：先验 p(θ)，红色：似然 p(y|θ)，绿色：后验 p(θ|y)。白色小圆为 MH 链在后验上随机游走的当前位置。

先验 / 似然 / 后验 PDF

MCMC trace plot — θ vs iteration

理论与主要公式

$$p(\theta|y) \propto p(y|\theta)\,p(\theta),\qquad \alpha = \min\left(1,\, \frac{\pi(\theta^{*})}{\pi(\theta_t)}\right)$$

π 为未归一化的后验密度。提案 θ* 以概率 α 被接受，足够长的链给出来自 π 的样本（Metropolis–Hastings）。

$$\frac{1}{\sigma_{\text{post}}^{2}} = \frac{1}{\sigma_p^{2}} + \frac{N}{\sigma_y^{2}},\qquad \mu_{\text{post}} = \sigma_{\text{post}}^{2}\!\left(\frac{\mu_p}{\sigma_p^{2}} + \frac{N\,\bar{y}}{\sigma_y^{2}}\right)$$

正态-正态共轭模型的解析后验。N 增大时 1/σ_post² 线性增长，σ_post 按 1/√N 收缩。

$$\text{CI}_{95} = \mu_{\text{post}} \pm 1.96\,\sigma_{\text{post}},\qquad I = \log\!\frac{\sigma_p}{\sigma_{\text{post}}}\ [\text{nats}]$$

95% 可信区间与信息增益 I（先验到后验中减少的不确定性，单位 nats）。

贝叶斯校准（MCMC）— Metropolis–Hastings

🙋

贝叶斯校准就是把 CAE 模型对齐到实验数据吧？跟常规的最小二乘有什么不同？

🎓

问得好。最小二乘（最大似然）只返回"最贴合数据的那一个 θ"，而贝叶斯把 θ 当作一个概率分布来对待。先验 p(θ) 注入过去知识或合理范围，再乘以观测的似然 p(y|θ)，得到后验 p(θ|y)。结果不是"一个点"而是"带离散的分布"，因此可以自然地谈设计裕度与置信区间。例如把材料的杨氏模量估计为 200 ± 5 GPa，后续应力计算就能把这个带状信息原汁原味地传下去。

🙋

明白了。那为什么还要 MCMC 出场？后验分布不能直接写成解析式吗？

🎓

像本工具这种"共轭"情形（正态先验 + 正态似然）后验确实有闭式，所以解析后验（绿色曲线）和 MCMC 估计能完美吻合。但真实的 CAE 模型——比如 k-ε 湍流常数、非线性硬化律——评估 p(y|θ) 需要跑一次完整数值仿真，闭式根本写不出来。Metropolis–Hastings 只需要 π(θ*)/π(θ_t) 的比值，对任何黑盒似然都能跑，所以成了万能工具。

🙋

移动 MH 提案宽度 step 时，接受率变化剧烈。最优值怎么定？

🎓

经典结果是 Roberts–Gelman–Gilks (1997)：一维约 44%、高维约 23.4% 最优。经验法则是 σ_proposal ≈ 2.38·σ_post。step 过大，接受率掉到接近 0%，链就冻住；过小则接受率高但样本自相关强，ESS 反而下降。瞄准 30〜50% 区间调参就是实务感觉。默认 step=0.15 对应 σ_post ≈ 0.07，比值约 2.1，正好落在最优点附近。

🙋

最后一个：信息增益是什么？以 nats 为单位还是第一次见。

🎓

它用信息量单位测量先验与后验之间的差距。这里用近似 log(σ_p/σ_post)，表示"看到数据后不确定性缩小了多少"。底取 e 得 nats（底取 2 则得 bits）。默认参数下 log(2/0.0707) ≈ 3.34 nats ≈ 4.8 bits，相当于"实验给出大约 5 比特关于 θ 的新信息"。比较不同校准实验的"信息价值"时非常方便。

常见问题

最小二乘（最大似然）只给出参数的一个最优点估计，而贝叶斯校准把参数 θ 看作一个概率分布：把表达先验知识的 p(θ) 与表达观测的 p(y|θ) 相乘，得到后验分布 p(θ|y)。这样不仅有最优值，还能定量描述其离散程度以及"取其他值的概率"。CAE 校准常面对数据少、噪声大的情形，单点估计容易过度自信。在本工具中将观测数 N 从 5 滑到 500，可以直接看到后验标准差 σ_post 按 1/√N 的规律收缩。

根据 Roberts-Gelman-Gilks (1997) 的渐近最优化结果，一维下最优接受率约为 44%，高维下约为 23.4%。经验法则是把提案标准差取为后验标准差的约 2.38 倍。本工具用 bell-shape 模拟了这一规律：step 过大时接受率几乎降为 0%，链子卡住不动；step 过小时接受率虽高但自相关强，有效样本数 ESS 会下降。实务做法是反复调节 step，使接受率落在 30〜50% 之间。

MCMC 的初值往往远离后验高密度区，最初的几百到几千步是"向高密度区漂移"的过渡样本。如果把它们混入统计，会让后验均值与方差产生偏差，因此需要丢弃，这就是 burn-in。本工具默认 500 步对于弱信息先验与数千样本的链子是较保守的取值。观察 trace plot，确认左端的漂移已消失即可。严谨场合需要并行多条链并用 Gelman-Rubin 的 R̂ < 1.01 判断收敛。

MCMC 样本之间存在自相关，因此 5000 个原始样本"折算为独立样本数"可能只有几百甚至几十，这就是 ESS。后验均值的标准误为 σ_post/√ESS，当 ESS 跌破 100 时，估计值的标准误就超过后验标准差的 10%，结果的小数点第二位以下就不能信任了。提高 ESS 的三种途径是：(1) 优化提案宽度；(2) 增加样本数；(3) 切换到 Hamiltonian Monte Carlo / NUTS 等更高效的采样器。

实际应用

CFD／FEM 代码的参数校准：湍流模型的常数（k-ε 的 C_μ、Cε1、Cε2 等）、非线性材料的硬化律参数、混凝土损伤模型的极限应变等，那些过去由分析师"拍脑袋"决定的常数，可以用实验数据通过后验分布来估计。NASA 与美国能源部已把基于 MCMC 的贝叶斯校准列为 VVUQ（Verification, Validation, Uncertainty Quantification）流程的必备步骤。本工具的正态-正态模型是其中最简单的共轭入门情形。

气候与地球科学模型校准：气候敏感度、云反馈系数等无法直接观测的参数，常用 MCMC 通过历史气温趋势或卫星观测来估计其后验。Hadley Centre 与 NOAA 的评估报告中后验分布的直方图是标配。在观测稀少、先验影响很大的领域，本工具中滑动 σ_p 可以直观看到先验对后验的拉动力度。

贝叶斯机器学习：贝叶斯神经网络（BNN）、高斯过程的超参数推断、A/B 测试中效应量估计等都依赖 MCMC。现在的主流是 HMC／NUTS（Stan、PyMC、NumPyro），但其底层的接受／拒绝思想正是本工具所演示的 Metropolis–Hastings。监控接受率与 ESS——"链子的健康检查"——在 HMC 中同样适用。

药代动力学（PK/PD）建模：用少量受试者数据为体内动态的房室模型估计带个体差异的后验，这一做法已成为行业标准（分层贝叶斯 + MCMC）。FDA 审查资料中后验可信区间是必填项。本工具中"观测数 N 越小、σ_post 越大"的现象，恰好对应受试者数量少的 I 期临床试验的结构。

常见误解与注意点

最大的坑是把"接受率高"等同于"采样器好"。把提案宽度 step 调得很小，接受率轻松超过 90%，但每步样本几乎不动，5000 个样本折算成独立样本可能只有几十个。在本工具中试试 step=0.05，看 ESS 怎样塌掉。接受率必须与 ESS 一起看，绝不能单独评估。

其次是"宽且弱信息的先验就一定客观"的误区。先验过宽时，少量数据会把后验"拉得很厉害"，离群值与极端值出现的频率会上升；反之，强信息的先验（σ_p 取小）表达对模型的自信，但也会掩盖模型与数据的不一致。在本工具中把 N=5、σ_p=0.5 与 N=5、σ_p=10 对比，可以亲身感受先验的影响。专业做法是用多个先验做敏感性分析，而不是一锤定音。

最后是只凭一条链的外观判断 MCMC 收敛。一条看似"稳定"的链可能只是停在了某个模式上；换个初值跑第二条链，可能收敛到完全不同的区域（典型于多峰后验）。生产场景下要并行多条链，要求每个参数的 Gelman–Rubin R̂ < 1.01。本工具是一维一链的概念学习器，只显示 R̂ 代理值。在真实项目中请务必用 PyMC 或 Stan 的 summary() 检查 R̂、ESS 与 MCSE。

贝叶斯校准（MCMC）— Metropolis–Hastings

常见问题

实际应用

常见误解与注意点

使用指南

具体计算示例

实务注意事项