贝叶斯校准
理论与物理
什么是贝叶斯校准
老师，贝叶斯校准是用实验数据来估计模型的参数吗？和普通的优化有什么不同呢？
问得好。简单来说，贝叶斯校准的特点是将参数估计为“概率分布”而非“单一数值”。通常的最小二乘法或优化是进行点估计，给出“最符合实验数据的值就是这个！”，而贝叶斯方法则返回一个分布形式的答案，告诉我们“参数以这样的概率落在这个范围内”。
返回分布形式的答案有什么好处呢？
举个例子，假设在汽车碰撞分析中估计FEM模型的杨氏模量。优化方法只会给出 $E = 210\,\text{GPa}$ 这一个值。但贝叶斯方法会给出“$E$ 是均值为 $210\,\text{GPa}$、标准差为 $5\,\text{GPa}$ 的分布”。这种不确定性信息可以直接传递给后续的可靠性分析。
更重要的是，它能够正式地融入工程师的经验知识（先验分布）。可以结合“这种材料的杨氏模量大概在 $200 \sim 220\,\text{GPa}$ 左右”这样的先验知识和实验数据，共同来缩小参数范围。这是贝叶斯方法最大的优势。
能把经验知识融入数学公式，这太有实用价值了！但这具体是怎样的数学呢？
贝叶斯定理与基本方程
贝叶斯校准的数学骨架是贝叶斯定理。对于参数向量 $\boldsymbol{\theta}$ 和观测数据 $\mathbf{D}$：
各项的物理意义
- $p(\boldsymbol{\theta}|\mathbf{D})$（后验分布）：看到实验数据后参数的概率分布。这是校准的最终输出。
- $p(\mathbf{D}|\boldsymbol{\theta})$（似然函数）：在给定参数值下获得实验数据的概率。衡量模型预测与实验数据的“偏离程度”。
- $p(\boldsymbol{\theta})$（先验分布）：在看到数据前关于参数的知识。反映文献值、经验、物理约束。
- $p(\mathbf{D})$（证据/边缘似然）：归一化常数。用于模型比较，但在参数估计中可作为常数忽略。
似然函数具体是怎么写的呢？
考虑最简单的情况。假设对于 $N$ 个实验数据 $d_i$，其与模型预测 $M(\mathbf{x}_i, \boldsymbol{\theta})$ 的差值服从独立同分布的正态误差：
这里 $\sigma$ 是观测噪声的标准差。可以将其作为已知固定值，也可以和 $\boldsymbol{\theta}$ 一起用贝叶斯方法估计。在实际工作中，将 $\sigma$ 也作为估计对象通常更稳健，因为很多时候无法准确知道误差的大小。
但是老师，FEM模型在一定程度上近似了物理吧？假设模型与实验的差异仅用“正态误差”就能解释，这个假设是否过于宽松了？
问得很尖锐！这正是Kennedy-O'Hagan框架所要系统处理的问题。
Kennedy-O'Hagan框架
Kennedy和O'Hagan于2001年提出的框架，已成为CAE贝叶斯校准事实上的标准。其核心是将实验值分解如下：
Kennedy-O'Hagan各项的含义
- $M(\mathbf{x}, \boldsymbol{\theta})$（仿真输出）：给定输入条件 $\mathbf{x}$ 和校准参数 $\boldsymbol{\theta}$ 时，FEM/CFD等的计算结果。
- $\delta(\mathbf{x})$（模型失配项）：无论怎样调整参数都无法弥补的模型本身的结构性误差。源于网格离散化误差、被忽略的物理现象、理想化的边界条件等。通常用高斯过程（GP）表示。
- $\varepsilon$（观测误差）：实验测量伴随的随机噪声。通常假设 $\varepsilon \sim \mathcal{N}(0, \sigma_\varepsilon^2)$。
加入 $\delta(\mathbf{x})$ 和不加入，具体会有什么不同呢？
如果不加入 $\delta$，模型的结构性误差就会被全部推给参数 $\boldsymbol{\theta}$ 来承担。例如，实际上是由于边界条件不完善导致模型与实验不符，却通过将杨氏模量调整到极端值来强行匹配。这就是所谓的参数补偿问题，会导致估计出物理上无意义的参数值。
通过显式地引入 $\delta(\mathbf{x})$，可以将“模型无法再现的部分”分离出来，从而将参数估计保持在物理上合理的范围内。
原来如此…！这是在承认模型局限性的基础上进行估计啊。但是 $\delta$ 本身是如何确定的呢？
$\delta(\mathbf{x})$ 是作为高斯过程（GP）从数据中学习的。也就是说，$\delta$ 的形状不是事先决定的，而是通过贝叶斯推理从实验数据与模型预测的残差模式中自动估计出来的。GP的核函数（例如平方指数核）的超参数也会一同被估计。
先验分布的选择方法
先验分布是怎么决定的呢？会不会变得主观和随意？
先验分布的设置确实是贝叶斯方法中最有争议的一点。在实际工作中，通常根据以下情况区分使用：
- 信息性先验分布：基于材料数据库的值、过去的实验结果、手册值来设定。例：结构钢的杨氏模量设为 $E \sim \mathcal{N}(210, 5^2)\,[\text{GPa}]$。
- 弱信息性先验分布：宽松地约束物理上合理的范围。例：摩擦系数设为 $\mu \sim \text{Uniform}(0.05, 0.8)$。
- 无信息性先验分布：如Jeffreys先验等，希望让数据说话时使用。但在高维情况下后验分布收敛会变慢。
在实际现场，弱信息性先验分布用得比较多。这是一种“排除物理上不可能的领域，但让数据有足够影响力”的平衡。
先验分布的选择有时会极大地改变结果吗？
当数据量少时，先验分布的影响会很大。反之，如果有足够的数据，即使先验分布有些不同，后验分布也会收敛到几乎相同的地方。这称为贝叶斯的渐近一致性。因此，在实际操作中，必须进行先验分布的敏感性分析（改变先验分布，确认结果在多大程度上变化）。
数值解法与实现
MCMC法基础
贝叶斯定理本身很简单，但实际计算后验分布似乎很困难吧？参数多的话积分就做不了了吧？
说得对。后验分布的归一化常数 $p(\mathbf{D})$ 是高维积分，无法解析求解。因此要使用马尔可夫链蒙特卡洛（MCMC）法。这是一种直接从后验分布生成样本的方法，增加样本数就能以任意精度近似后验分布。
MCMC是什么机制呢？
最基本的Metropolis-Hastings算法流程如下：
- 设定初始值 $\boldsymbol{\theta}^{(0)}$
- 从提议分布 $q(\boldsymbol{\theta}^*|\boldsymbol{\theta}^{(t)})$ 生成候选点 $\boldsymbol{\theta}^*$
- 计算接受概率：
以概率 $\alpha$ 接受候选 $\boldsymbol{\theta}^*$，以 $(1-\alpha)$ 拒绝并停留在当前值。将此过程重复数千至数万次，样本序列就会收敛到后验分布。
可以想象成在山中随机游走，并在海拔高（似然×先验概率大）的地方停留更久的机制。
哦，原来如此！但是FEM计算一次要几分钟到几小时吧？运行几万次现实吗？
问得很尖锐。这正是CAE中使用贝叶斯校准的最大瓶颈。对策是稍后要详细讲的代理模型，不过我们先整理一下各MCMC算法的特点。
MCMC算法比较
| 算法 | 特点 | 适用场景 | 典型接受率 |
|---|---|---|---|
| Metropolis-Hastings | 最基本。实现简单 | 低~中维度（~10个参数） | 20〜50% |
| Adaptive Metropolis (AM) | 自动调整协方差矩阵 | 中维度。无需调参 | 23%（最优） |
| DRAM | 延迟拒绝 + AM。拒绝时重新提议 | 中维度。对多峰性有一定应对能力 | 30〜50% |
| 哈密顿蒙特卡洛 (HMC) | 利用梯度信息。高维也高效 | 高维度（100+）。可微分模型 | 65〜90% |
| NUTS | HMC的自动调参版。Stan标准 | 高维度。无需手动干预 | 80〜95% |
| Ensemble Sampler (emcee) | 多个行走器并行探索 | 中维度。易于并行化 | 20〜50% |
| TMCMC / SMC | 序贯蒙特卡洛。对多峰性强 | 多峰后验分布。模型比较 | — |
有这么多啊！CAE中常用的是哪种呢？
Sandia国家实验室的Dakota中DRAM是标准实现，在实际工作中最常用。如果能获取梯度信息，HMC/NUTS的效率则极高。用Python轻松实现的话，推荐PyMC（基于NUTS）或emcee。
收敛诊断
怎么判断MCMC“收敛”了呢？
这是非常重要的一点。如果使用了未收敛的MCMC结果，会得到完全没有意义的后验分布。标准的收敛诊断如下：
- $\hat{R}$（Gelman-Rubin统计量）：多个独立链的方差比。若 $\hat{R} < 1.01$（严格来说是 $< 1.05$）则判断为收敛。
- 迹图：绘制样本序列的时间序列图。理想的混合良好的模式是“毛毛虫状”。若有趋势或跳跃则未收敛。
- 有效样本大小 (ESS)：考虑自相关后的实质样本数。目标为 ESS $> 400$（每个参数）。
- 老化期剔除：舍弃初期瞬态样本（通常为链长的10〜50%）。
确认 $\hat{R} < 1.01$，迹图稳定，且ESS足够的话，就OK了是吧。
基于代理模型的加速
刚才说到“FEM计算运行几万次很困难”，那实际上是怎么解决的呢？
最通用的解决方案是代理模型（替代模型）。用快速的近似模型替代FEM的输入输出关系。主要方法有：
- 高斯过程回归（GP/Kriging）：与Kennedy-O'Hagan框架亲和性高。也能给出预测的不确定性。
- 多项式混沌展开（PCE）：用正交多项式展开输入参数的概率分布。可同时用于不确定性传播。
- 神经网络：在高维、强非线性情况下有效。但需要数据量。
例如，执行FEM 100〜500次来学习GP，那么在MCMC的每一步就只需要GP预测（毫秒级）。理论上可以将FEM的计算次数从 $10^4$ 减少到 $10^2$。
如果代理模型的精度不好，校准结果也会出错吧？
没错。所以必须验证代理模型。具体来说，通过留一法交叉验证来确认预测精度。目标是 $Q^2 > 0.99$。另一个策略是序贯代理模型更新。在MCMC过程中，在后验分布的高概率区域添加FEM计算点，逐步精细化代理模型。这样可以兼顾效率和精度。
实践指南
校准工作流程
请告诉我实际进行贝叶斯校准时的步骤！
なった
詳しく
報告