带有似然温度调节和并行实现的序贯蒙特卡洛方法在不确定性量化中的应用

《AIChE Journal?AIChE》：Sequential Monte Carlo with likelihood tempering and parallel implementation for uncertainty quantification

　　摘要

贝叶斯估计能够实现不确定性的量化，但其解析实现往往难以处理。作为一种近似方法，马尔可夫链蒙特卡洛（MCMC）方法被广泛使用，尽管由于需要频繁评估似然函数而带来了较高的计算成本。作为另一种方法，序贯蒙特卡洛（SMC）方法在时间序列数据方面表现良好，而其在非时间序列数据上的优势（通过似然温度调节）尚未得到充分探索。在本研究中，我们实现了一个带有似然温度调节的SMC算法，并通过两种类型的问题将其与MCMC进行了比较：玩具问题和使用真实实验数据的甲烷化反应器模型的参数不确定性量化问题。此外，我们还评估了并行化对SMC性能的影响。结果表明，对于参数相关性强且需要计算成本高昂的模拟的似然评估问题，带有似然温度调节的SMC在计算效率和稳定性方面表现出显著优势。

1 引言

数学模型被广泛用于化学过程设计和天气预报等目的。它们的参数对结果有很大影响，但这些参数通常是未知的或仅能近似知道，这限制了模拟的准确性。因此，需要进行参数估计以获得可靠的参数并建立适合预测应用的数学模型。存在多种参数估计方法。其中最简单的方法是最小二乘法。然而，这种方法仅提供点估计。基于单一参数集的模型只能产生一个预测结果，无法评估模型的可靠性。为了构建稳健的数学模型，需要能够量化不确定性的参数估计方法。一种解决方案是基于贝叶斯定理的估计，即贝叶斯估计，其中参数的先验分布通过观测数据更新以获得后验分布。通过将先验信息与观测数据结合到估计过程中，可以有效地获得参数的不确定性。然而，解析执行贝叶斯估计通常是不可行的或不可能的。为了解决这个问题，使用了近似估计方法。一种代表性的方法是马尔可夫链蒙特卡洛（MCMC）方法。MCMC方法是一种利用马尔可夫链的基本属性的方法，同时确保收敛性、非周期性和详细平衡。在这种方法中，样本沿着马尔可夫链依次生成，每个样本对应于数学模型中的一组参数值。从初始参数值开始，算法提出新的候选值，并根据指定的标准接受或拒绝它们，从而通过采样近似目标概率分布。许多研究已成功将MCMC应用于各种领域的数学模型中的参数估计。尽管如此，MCMC方法的一个主要缺点是计算成本高。作为替代方案，序贯蒙特卡洛（SMC）方法提供了类似于MCMC的近似贝叶斯推断。SMC方法最初是为状态空间模型中的序贯状态估计开发的，现在被广泛用于时间序列预测。SMC方法包括四个主要步骤：初始化、更新、权重计算和重采样。MCMC和SMC之间的关键区别在于目标分布的近似方式。MCMC通过沿着马尔可夫链的顺序采样来构建后验近似，而SMC使用基于总体的目标分布表示，并在每次迭代时更新和重新加权。SMC和MCMC之间最显著的区别在于并行性。对于这两种方法，限制计算速度的步骤是似然计算，它占据了总体计算成本的绝大部分。在MCMC中，模拟和似然评估是在单个马尔可夫链内依次进行的，使得在这个计算成本高昂的步骤中无法实现并行化。虽然可以并行运行多个独立的马尔可夫链，但链级并行性受到每个链的收敛性评估的实际考虑的影响。相比之下，SMC可以同时为多个粒子执行似然计算，从而实现并行化，从而有可能更快地完成估计。随着并行计算的不断进步，预计SMC在计算效率方面的优势将会增加。尽管有这些优势，但使用SMC进行参数估计的报告仍然不如基于MCMC的方法多。SMC包括两种主要变体：数据温度调节和似然温度调节。数据温度调节逐渐引入数据以转换分布，使其非常适合时间序列应用，尽管估计精度取决于数据包含的顺序。相比之下，似然温度调节（也称为桥接方法）逐渐增加似然的影响，以实现从先验分布到后验分布的平滑过渡，并且主要应用于非时间序列数据的研究中。然而，利用似然温度调节的研究相对较少，对其性能的定量评估尤其稀缺。因此，系统地评估带有似然温度调节的SMC的有效性和特性具有学术意义。在本研究中，我们实现了一个带有似然温度调节的SMC算法，这在之前的工作中很少被研究，并对其实用性和性能进行了定量评估。通过将其与代表性的贝叶斯估计技术MCMC进行比较，我们展示了SMC在不确定性量化精度和计算效率方面的优势。我们还展示了并行化在降低SMC计算成本方面的有效性。最后，我们将所提出的方法不仅应用于玩具问题，还应用于实际化学工程系统——甲烷化过程模型的参数估计和不确定性量化，并证明了SMC在实际应用中的有效性。

2 背景和方法

2.1 贝叶斯推断

在本研究中，贝叶斯估计用于参数估计。贝叶斯估计是一种基于贝叶斯定理的方法，可以定量估计观测结果背后的参数概率分布。使用贝叶斯估计，参数可以作为概率分布而不是单一值来估计。此外，在估计开始时指定的先验分布允许将先验知识和经验信息纳入分析。贝叶斯定理可以用方程（1）表示。

（1）

其中 是参数， 是数据； 是先验分布； 是似然； 是后验分布。此外， 表示边际似然。在本研究中，实验数据 和模型近似 如方程（2）所示相关。

（2）

其中 表示假设为独立同分布的观测误差，遵循正态分布：

（3）

这里， 是与实验误差相关的标准差； 是数据点的数量； 假设独立遵循相同的正态分布。每个实验条件的似然可以使用方程（4）计算。

（4）

通过将每个实验条件的似然相乘，似然可以用方程（5）表示。

（5）

在化学过程中估计参数时，评估 通常计算成本高昂，导致评估方程（5）的计算成本很高。在某些示例中，这种计算成本尤其显著，包括分布式参数系统或动态过程，例如本研究中讨论的那些。因此，减少 的评估次数非常重要。在本研究中，测量噪声被建模为独立的和高斯的（方程（4））。对于这里考虑的甲烷化反应器，这种建模选择是合理的，因为观测是在稳态操作条件下获得的，每个数据点对应于一个独立进行的实验。此外，测量误差来自多个来源，如流量测量、温度控制和气体成分分析，其综合效应可以根据中心极限定理合理地用高斯分布近似。在化学反应器模型的参数估计研究中，也广泛采用了类似的假设，特别是当在稳态条件下使用最小二乘或基于似然的公式时，提供了一种实用且易于处理的似然公式。

2.2 SMC方法

由于解析执行贝叶斯估计可能具有挑战性，因此开发了几种近似后验分布的方法。SMC方法最初是为状态空间模型中的状态变量序贯估计开发的，现在被广泛用于时间序列数据的预测。SMC方法通过使用多个粒子跟踪变量来估计后验分布，这些粒子基于从观测数据获得的信息。SMC方法的一个显著特点是易于并行化。与稍后描述的MCMC方法不同，解决每个粒子的模型不会影响其他粒子。因此，每个粒子可以独立处理。通过将粒子分配给多个核心进行处理，可以轻松实现并行化。此外，可以根据它们的先验分布在整个参数空间中对粒子进行采样，并根据每个点的似然评估进行平衡，从而实现稳定采样，而不会陷入参数空间的特定区域。这种行为在第3节中进行了演示。在SMC方法中，使用两种方法构建中间分布，这是通过修改贝叶斯定理（方程（1）中的似然项来实现的，以便逐渐引入观测数据的贡献。这种构建使得从先验分布到后验分布的过渡平滑，并构成了SMC中顺序更新过程的基础。数据温度调节逐渐将数据添加到似然计算中以更新中间分布。另一方面，似然温度调节逐渐增加似然对中间分布的影响。数据温度调节和似然温度调节之间的区别在于中间分布的计算方式。在数据温度调节中，数据集被分成子集并依次纳入似然函数中。然后得到中间分布：

（6）

这里， 表示完整的观测数据集；索引k表示单个观测； 表示第k个数据点。项 表示给定参数向量 的第k个观测的似然贡献。在温度调节阶段， 表示纳入后验分布的观测数量，其中， 表示可用观测的总数。因此，中间分布 对应于基于前 k 个观测的后验分布。在数据温度调节中，每个中间分布 可以解释为基于观测数据子集的后验分布。相比之下，在似然温度调节中，中间分布的获得方式如下：

（7）

这里， 表示估计步骤的数量， 是估计过程终止所需的估计步骤总数。此外， 称为温度调节因子，随着估计的进行而单调增加从0到1，当 时，它减少到方程（1）中的分子。因此， 表示第i步的中间分布。应该注意的是，与数据温度调节不同，似然温度调节不将中间分布视为后验分布。图1展示了带有似然温度调节的SMC算法的概述，它包括四个主要步骤。在步骤I（初始化）中，定义了一个先验分布，并从这个分布中抽取了一组初始样本。接下来，在步骤II（似然计算）中，评估每个样本的似然，并根据每个样本解释观测数据的程度分配权重，使得似然值较高的样本获得较大的权重。在步骤III（重采样）中，具有相对较大权重的样本被概率性地复制，而权重较小的样本更有可能被丢弃。最后，在步骤IV（变异）中，每个样本被分配一个来自提议分布的提议参数值，并进行随机更新，完成一次迭代后似然计算阶段再次开始。

图1 显示了带有似然温度调节的序贯蒙特卡洛（SMC）算法的示意图。左侧面板显示了算法的工作流程，而右侧面板展示了粒子在温度调节阶段的相应演变。每行水平线代表在给定温度调节阶段获得的样本分布，对应于温度调节参数的特定值，而水平轴代表参数空间θ。完成这一系列步骤后，温度调节因子 增加，算法返回到步骤II。通过逐渐增加 ，构建了一系列连续桥接先验和后验分布的中间分布。当突变步骤完成时，估计过程终止，此时收集的样本集合提供了目标后验分布的近似值。图1的左侧展示了SMC算法的流程，而右侧则示意性地描绘了样本在连续的退火阶段中的演变过程。随着算法通过似然计算、重采样和突变的进行，样本逐渐集中在后验概率较高的区域，直观地展示了后验分布是如何从先验分布构建而来的。图2概念性地比较了SMC方法和马尔可夫链蒙特卡洛（MCMC）方法在处理复杂后验分布时的采样能力，其中图2A展示了典型的单链MCMC过程，图2B展示了SMC方法。在典型的单链MCMC设置中，样本是沿着马尔可夫链顺序生成的，对后验分布的探索可能受到初始状态和随机种子的影响，从而导致对参数空间的探索效率低下。相比之下，SMC方法维护了一个初始时覆盖整个参数空间的样本群体，并通过退火、重采样和突变逐渐集中在高概率区域，遵循图1中所示的迭代过程。因此，SMC能够更稳健地探索复杂的后验分布。由于在似然评估和突变步骤中每个粒子都是独立处理的，SMC也允许直接进行并行化。该图旨在说明采样行为上的概念性差异，而不是提供定量的性能比较。

图2展示了在二维参数空间中（A）MCMC和（B）SMC的采样方法。在MCMC示例中，黑色圆圈和正方形代表来自不同初始值的样本，箭头表示沿马尔可夫链的转移；而在SMC示例中，灰色圆圈代表从先验分布中抽取的粒子，黑色圆圈代表近似后验分布的粒子。

2.3 SMC算法
整体SMC算法在算法1中呈现，该算法遵循图1中所示的概念流程图。重采样步骤在算法2（附录S1的第S2节）中有详细说明。该算法的源代码可以在我们的GitHub仓库中找到。

算法1. 带有似然退火的SMC方法

在下文中，我们将介绍算法1中出现的符号和超参数。第m个粒子的参数向量表示为，其中每个粒子代表参数空间中的一个候选值。这些粒子通过它们在参数空间中的分布共同近似后验分布。基于粒子集合，协方差矩阵用于描述后验的分布和关联结构，并按以下公式计算：

这个协方差矩阵在算法1（第24行和第29行）中被使用。在算法1（第37行）中，我们使用了单位阶跃函数。单位阶跃函数表示为

算法中使用的超参数定义如下：表示粒子的总数；是退火因子的最大增量；是动态调整增量衰减率。是突变步骤中用于评估提议移动效率的接受率阈值。在算法1中，似然评估步骤在多个计算核心上并行进行，其中表示可用核心的数量；是标识每个核心的索引；是分配给第c个核心的粒子集合。这种分配可以在执行过程中动态确定，以实现多个核心之间的负载平衡，并在似然评估步骤（第12行和第31行）中使用。在本研究中开发的SMC方法中，退火因子逐渐增加到1，从而逐步将数据中的信息引入中间分布。这里，增量的选择必须谨慎。如果信息反映得太强烈，就会发生一种称为退化现象——即重采样后，大多数粒子的权重可以忽略不计，只有少数粒子主导近似结果，导致多样性差和后验估计不稳定。另一方面，如果增量太小，则收敛所需的时间会变长。为了找到合适的增量，我们检查了有效样本大小（ESS），如算法1（第17行）中所评估的。

2.4 与MCMC方法的比较
MCMC方法是一种利用马尔可夫链的性质来生成逐渐遵循给定后验分布的样本序列的技术。MCMC方法和SMC方法之间的主要区别在于，MCMC方法是顺序生成样本的。由于这一特性，MCMC算法在单个马尔可夫链内的并行化本质上是困难的，因此并行执行主要限于链级并行性或特定组件（如下所述，如似然计算）。MCMC的基本算法包括Metropolis–Hastings（MH）方法和切片采样。这些算法顺序生成候选参数值，并根据它们的接受概率确定下一个样本。一旦马尔可夫链通过足够的迭代达到平衡，生成的样本就被用来近似后验分布。在本研究中，我们使用Python库PyMC中的这两种算法与SMC进行比较。

我们比较了PyMC中实现的两种无梯度MCMC采样器：Slice采样器和MH采样器。Slice采样器使用其默认配置，该配置会自适应地调整步长，无需手动调整。对于MH采样，我们使用了PyMC中的默认Metropolis实现，基于高斯随机游走提议和默认的缩放参数。这两种采样器仅依赖于似然评估，不需要梯度信息，确保了与本研究中考虑的仅基于似然的设置一致。除了第3.2.1节中描述的实验外，MCMC采样器的燃烧期（热身期）设置为PyMC提供的默认值；在第3.2.1节中，燃烧长度是明确指定的，以适应实验设置。尽管基于梯度的采样器（如No-U-Turn Sampler (NUTS)）广泛用于连续和相关参数空间，但它们的适用性严重依赖于对数后验密度的准确梯度的可用性。在本研究中，似然评估是通过外部数值ODE/DAE求解器进行的，该求解器被视为黑盒模拟器，自动微分并不容易获得。实现可靠的敏感性信息（例如，伴随或前向敏感性）需要超出本工作范围的大量方法论开发，而有限差分近似会显著增加计算成本并引入数值噪声，这已知会对哈密顿蒙特卡洛方法的性能和数值稳定性产生不利影响。因此，为了在相同的信息约束下进行公平比较，我们将MCMC比较限制在仅基于似然评估的无梯度采样器上，这与所提出的SMC框架的假设一致。在本研究中，我们关注MCMC中两种常见的并行化形式：链级并行性（多个独立的马尔可夫链同时执行）和似然级并行性（在数据点或实验条件上并行评估似然函数）。链级并行性通常用于诊断目的，而不是为了加速采样效率。当多个链并行运行时，每个链仍然必须顺序演化并生成足够长的轨迹以充分探索目标后验分布。收敛性通常使用Gelman–Rubin诊断和ESS等诊断指标进行评估，这两个指标都是在链级别定义的或在链之间聚合的。在MCMC中，ESS通过考虑每个马尔可夫链内的自相关性来量化等效的独立样本数量。重要的是，要获得接近1的可靠ESS值和足够大的ESS，通常需要每个链包含大量的燃烧后样本。因此，在固定的计算预算下，将样本分布在许多短链上可能导致链内混合不良和链级ESS较低，从而限制了通过链级并行性实现的实际性能提升。因此，链级并行性不能被视为与SMC中的基于粒子的并行性完全等效，在SMC中，使用粒子群体共同构建单一的后验近似。在本研究中，我们还通过运行多个独立的MCMC链（总样本数固定）在案例1.2中检验了链级并行性的实际效果。相比之下，似然级并行性利用了似然函数本身的结构，通过分解方程（5）中的求和来并行评估似然函数。这种并行形式与链级并行性互补，特别是在似然评估涉及模拟时特别有效，因为模拟占据了计算成本的主体。在本研究中，似然级并行性在案例2.1和案例2.2中用于MCMC，其中对每个实验条件独立评估似然。这种方法允许高效利用多个计算核心，而不改变MCMC算法的基本顺序结构。另一方面，SMC方法使用多个样本（粒子）来近似中间和后验分布，这些样本可以轻松分配到多个处理核心（算法1中的第12行和第31行）。使用每个粒子的参数进行评估是相互独立的。因此，当粒子数量超过核心数量时，完成一个粒子计算的核可以立即重新分配给下一个核心。这样，核心可以连续处理任务，而无需等待，从而实现计算资源的有效和连续利用。

2.5 量化分布的近似程度
我们量化了真实分布和采样结果之间的差异，两者都是一维的：

这里，sup表示所有x上的最大值（上确界）。这个指标直接测量了两个累积分布之间的最大差异，无需进行密度估计，因此对于有限样本来说是稳健的。尽管它通常用作Kolmogorov–Smirnov（KS）检验中的检验统计量，但在这里仅作为近似精度的定量度量使用，不进行假设检验。应当注意的是，当样本表现出强自相关性时，不应在正式的假设检验框架内解释KS统计量。

2.6 案例研究
在本研究中，我们讨论了几个数学模型，以基准测试SMC和MCMC方法。

2.6.1 玩具问题
为了测量SMC和MCMC的采样精度，我们考虑了以下五个由方程（12-16）–（12-16）表示的二元函数。案例1.1：变量间相关性弱的二元正态分布。
案例1.2：变量间相关性强的二元正态分布。
案例1.3：香蕉形状的函数。
案例1.4：环形函数。
案例1.5：多模态函数。
在本研究中，可以直接从分布范围内进行采样。每种方法从真实分布中采样的能力通过方程（11）进行评估。

2.6.2 甲烷化过程模型
除了上述玩具问题外，我们还考虑了一个化学工程中的实际问题。我们基于贝叶斯估计量化了CO2甲烷化过程中的参数不确定性。该模型和数据来自Tsuboi等人之前的研究。过程流程图在附录S1中提供。实验的详细信息如下：使用的观测数据是每种物质的出口流量。输入变量包括反应器入口处的气体组成、反应器入口处的压力、流速、夹套温度、催化剂用量以及催化剂颗粒大小。输出是出口处各组分的摩尔分数。这些摩尔分数是通过两台气相色谱仪获得的。所使用的催化剂填充在反应器内部，并层叠在中心的烧结过滤器上。使用氩气来稀释进料气体并跟踪体积变化。二氧化碳的甲烷化反应通过方程（17）所示的反应进行。

正如Tsuboi等人所展示的，该模型是在以下假设下开发的：

- 忽略了催化剂与气相之间的温差。
- 忽略了反应器径向方向上的浓度和温度分布。仅考虑二氧化碳的甲烷化反应，忽略了其他副反应，如一氧化碳的甲烷化。
- 反应速率由Langmuir-Hinshelwood模型描述。
- 应用了理想气体状态方程。
- 可以忽略孔隙扩散的影响。

物料平衡方程由以下公式给出：

其中，是床层空隙率；是时间；是达到稳态的时间；是反应器的长度；是气相中组分的浓度；是表面速度；是轴向扩散系数；是组分的化学计量系数；是反应速率。方程左侧的项代表各组分的积累，而右侧的项分别代表由于对流、床层内的轴向扩散和反应引起的变化。总物料平衡方程由下面的公式给出：

通过将理想气体状态方程应用于此方程，总物料平衡方程也可以表示为下面的公式：

能量平衡方程可以由以下公式给出：

其中，是气体密度；是催化剂密度；是气体的比热容；是催化剂的比热容；是有效轴向热导率；是反应焓；是夹套温度；是总传热系数；是床层直径。在本研究中，总传热系数的值是根据之前的研究估算的。方程左侧代表气相与催化剂之间的温度变化。由于催化剂颗粒足够小，假设气相与催化剂之间的传热速率足够快，气相和颗粒的温度会迅速达到平衡。右侧代表对流、床层内的轴向热传导、反应热以及向夹套管的热量散发。这里，通过以下公式估算：

模型中的反应速率方程与之前的研究中的相同。该模型假设二氧化碳的甲烷化反应速率由Sabatier反应中的正向和反向反应速率方程之差给出，每个反应速率都通过带有经验确定的指数的改进型Langmuir-Hinshelwood机制进行建模。该模型没有考虑一氧化碳的存在，因为检测到的一氧化碳量仅约为3%，被认为是可以忽略的。反应速率模型由以下公式表示：

此外，和分别是二氧化碳和水分子的吸附平衡常数，由以下Van't Hoff方程给出：

进一步地，和分别是二氧化碳和水分子的吸附平衡常数，由以下Van't Hoff方程给出：

通过解这些方程，可以获得稳态出口浓度、出口温度和出口流速，这些值在标准条件下转换为出口处各组分的体积流量：

方程（26）中的60和60因子仅用于单位转换：60因子将时间从秒转换为分钟，而因子将体积从转换为。为了数值求解该模型，我们使用了Python包assimulo，它提供了微分-代数方程（DAE）求解器的接口，并通过重复执行时间积分直到收敛来获得稳态解。具体来说，DAE系统是使用SUNDIALS套件中的IDA求解器解决的。采用这种方法是因为模型被构建为一组微分-代数方程，由于数值病态性，不能通过直接将时间导数设为零来获得稳态解。对于模型的每次评估，都进行时间积分直到达到稳态。我们设置足够大的时间，以便模型达到稳态。对于模拟所需的催化剂属性，我们使用了之前研究中报告的值。在参数估计中，我们仅将数学模型的稳态解与在稳态条件下获得的实验数据进行了比较。本研究中估计的参数如下九个参数：

从之前的研究中估计的参数总结在表S1中。在本研究中，我们使用了之前研究中报告的不同条件下进行的64次实验获得的实验数据。其中，进料气体是H2、CO2和Ar混合物的实验称为正向反应实验，而进料气体是CH4、H2O和Ar混合物的实验称为反向反应实验。总共进行了29次正向反应实验和35次反向反应实验，每种反应使用了两种条件作为验证数据，如Tsuboi等人所述。在这些实验中，改变了以下条件：夹套温度；入口压力；催化剂直径；每种催化剂的孔隙率；以及反应器长度。在标准条件下，获得了每种组分在出口处的体积流量数据。参数估计是通过计算与从方程（26）获得的计算结果之间的似然值来进行的。具体的实验条件显示在附录S1中。对于甲烷化过程，进行了两个案例研究来评估SMC方法和MCMC方法的性能。第一个案例使用模型生成的带有噪声的模拟（人工）数据，其中噪声来自均值为0、标准差为5的正态分布，相当于参考尺度值（100）的5%。第二个案例使用真实的实验数据。每个案例的具体条件如下所述。在案例2.1中，通过聚类60个实验条件生成了模拟数据。每个变量都进行了标准化，然后使用K-means方法进行聚类，并选择了30个具有附加噪声的代表性条件进行模拟。模拟数据的参数使用之前研究中获得的值进行设置，如表S1所示，噪声的标准差为。在案例2.2中，使用了所有60个实验数据点。估计的参数是表S1中显示的八个反应速率方程参数以及数据的标准差。对于参数的先验分布，除了之外，设置了正态分布，以之前研究中的估计值为中心，标准差为50%。对于，设置了均值为10的正态分布，标准差为50%。如下所示，与先验分布相比，后验分布的不确定性显著降低，表明实验数据包含了关于的足够信息，我们选择的先验分布没有影响SMC获得的后验分布。案例1和案例2的设置和关键特征总结在表1中。

表1. 案例1和案例2的总结。

3 结果与讨论

3.1 玩具问题

我们使用玩具问题来评估我们算法中不确定性量化的准确性、计算速度和超参数的敏感性分析（图S6）。每个问题都是使用以下方法进行的。

3.1.1 分布近似准确性的评估

对于案例1，我们分别使用SMC方法和MCMC方法（Slice采样和MH方法）进行了100次估计，样本量从1000到11,000不等。对于所有算法，每次估计的随机种子值都不同。然后我们计算了近似准确性，定义在方程（11）中，以比较采样准确性。图3显示了基于使用SMC和MCMC方法获得的样本的θ1的边际概率密度函数的近似准确性。近似结果以箱线图的形式绘制，其中水平轴上的Np代表样本数量。位于95%区间外的异常值没有显示在图中，箱内的白点代表平均值。较小的值表示分布近似的量化准确性更好。

图3显示了玩具问题中每种方法的分布近似准确性比较。红色：SMC方法，蓝色：Slice采样，黑色：MH方法。箱线图中箱中心的标记代表平均值。（A）案例1.1，（B）案例1.2，（C）案例1.3，（D）案例1.4，（E）案例1.5。在图3A的案例1.1（弱相关正态分布）中，SMC方法和Slice采样的最大值和最小值之间的差异、四分位数范围和平均值对于所有采样数量都是相似的。相比之下，MH方法在最大值和最小值、平均值和四分位数范围之间的差异较大。此外，在图3B的案例1.2（强相关正态分布）中，SMC方法的四分位数范围和平均值最小，其次是Slice采样和MH。这表明SMC方法可以从强相关分布中进行准确和一致的采样。在图3C的案例1.3（香蕉形函数）和图3D的案例1.4（环形函数）中也可以找到类似的趋势。最后，在图3E的案例1.5（多模态函数）中，SMC方法的性能明显优于Slice采样和MH。从这个案例研究中，我们发现即使样本数量较少，SMC方法也能近似具有挑战性的分布，并且近似准确性是一致的。图4显示了案例1.4（环形函数）的样本，每种方法采样了1000个点。在这个图中，样本以核密度估计的形式显示，它在采样点周围放置一个平滑的分布（核），并将它们相加以获得平均密度。案例1中其他情况的核密度估计总结在图S2-S5中。

图4显示了案例1.4（环形函数）的样本，其中（A）MH方法，（B）Slice采样，（C）SMC方法和（D）真实分布的核密度估计结果。在图4A中，对于MH方法，分布被分割在左侧，并且看起来是不连续的，采样点在顶部和底部局部分布。这表明通过这种方法获得的样本缺乏覆盖整个区域的能力。在图4B中，对于Slice采样，与图4D中的真实分布相比，观察到了一定程度的再现性。然而，在等高线中观察到了不对称和局部不规则性，并且样本倾向于偏向某些区域。与MH和Slice采样相比，图4C中的SMC方法更接近真实分布，其等高线比(A)和(B)中的更平滑。SMC方法的更好采样效果通过这个具有挑战性的函数得到了证明。

3.1.2 计算速度的评估
我们使用案例1.2（具有强相关性的正态分布）作为基准来比较计算成本。所有计算都在配备了AMD Ryzen Threadripper PRO 5995WX CPU（64个物理核心，64个逻辑处理器，2.7 GHz基础时钟频率）和128 GB RAM的工作站上执行。为了观察计算时间的差异，我们在每种方法的似然性评估过程中引入了30毫秒的人为延迟，以模拟底层数学模型计算成本较高的情况。在这个实验中，SMC方法的样本数量被设置为1000个，SMC的超参数被固定为。对于MCMC方法，样本的总数被固定为1000个。报告的MCMC方法的计算时间不包括烧录（热身）阶段，仅包括烧录后的采样阶段。在这项研究中，我们报告的是墙钟计算时间，而不是每秒的ESS。这个选择反映了在SMC和MCMC中有效样本大小有不同的含义。在SMC中，ESS主要衡量粒子简并度，而在MCMC中，它反映了马尔可夫链的自相关程度。因此，基于ESS的直接效率比较可能会产生误导。在这项研究中，我们没有考虑MCMC方法中的链级并行性作为主要性能比较的一部分。这个选择是为了确保在固定的总采样预算下进行公平比较，关注单次估计运行的效率，而不是通过将样本分布在许多短链上来获得的吞吐量。当固定的总样本数量分配给许多独立链时，每个链只能获得有限的样本数量，这可能不足以实现可靠的收敛和准确的后验估计。对于像Slice采样这样的无梯度采样器，需要足够的链内样本来确保对目标分布的充分探索。为了检验在固定采样预算下链级并行性的影响，我们还使用25个并行链对案例研究1.2进行了额外的Slice采样实验，每个链仅生成40个样本。详细结果在附录S1的表S2中提供。在这种设置下，由于每个链内的样本数量有限，标准的收敛诊断指标（如链级ESS值）和得到的值无法可靠地解释。尽管如此，均匀较低的链内ESS值和值与单位的偏差表明在这种采样配置下无法确信收敛。这些观察结果表明，为每个链分配少量样本不足以进行可靠的后验推断，从而支持了我们不在主要比较中包括链级并行性的决定。在匹配的计算预算下，系统地、完全公平地比较SMC和MCMC之间的并行化策略仍然是未来工作的一个重要课题。相比之下，SMC框架中采用的并行化策略与MCMC中的链级并行化有根本不同。如算法1（第12行和第31行）所示，SMC在每次迭代中对一组共同的粒子并行化似然性评估。这种方法提高了计算效率，同时产生的结果与单核执行得到的结果相当。图5A展示了在单核条件下MH方法、Slice采样和SMC方法的平均计算时间比较。垂直轴上的平均时间代表使用不同随机种子进行的五次独立估计的平均计算时间，通过MH方法的平均计算时间进行了归一化。在这种归一化下，SMC方法的计算时间大约是MH方法的2.3倍，而Slice采样需要的时间大约是MH方法的4.3倍。

3.2 甲烷化过程模型
关于第2.6.2节中描述的甲烷化过程模型的参数估计，我们比较了SMC方法和MCMC方法（Slice采样）。由于即使经过数十小时的估计，分布估计似乎也没有收敛到后验分布，因此没有考虑MH方法。

3.2.1 使用模拟数据的计算效率比较
案例2.1比较了每种方法的估计精度和计算时间，以定量评估计算效率的差异。SMC方法的超参数选择如下：，和。对于Slice采样，前2500个样本被丢弃作为烧录期。图6使用小提琴图展示了使用SMC和MCMC方法获得的采样结果。对于SMC方法，采用了粒子级并行性，即在每次迭代中，对于一组共同的粒子，同时在多个处理器上执行似然性评估，如算法1（第12行和第31行）所示。为了检验独立运行之间的变异性，展示了三个使用不同随机种子（SMC1–3）的独立估计；改变随机种子会影响算法1中的随机操作，包括粒子初始化和变异步骤（第4行、第29行和第34行）。同样，对于MCMC方法，展示了三个独立的估计（MCMC1–3），每个估计包含从不同初始值和随机种子生成的5000个样本。对于MCMC，采用了似然级并行性，其中对应于不同实验条件的似然贡献被并行评估，而马尔可夫链本身是顺序演化的。在这两种方法中，计算都使用了相同数量的计算核心（30个核心），以确保在计算努力方面的公平比较。MAP估计是通过使用差分进化算法最小化负对数后验来获得的，这是一种基于种群的全局随机优化方法。

3.2.1 使用模拟数据的计算效率比较
案例2.1比较了每种方法的估计精度和计算时间，以定量评估计算效率的差异。SMC方法的超参数选择如下：，和。对于Slice采样，前2500个样本被丢弃作为烧录期。图6使用小提琴图展示了使用SMC和MCMC方法获得的采样结果。对于SMC方法，采用了粒子级并行性，即在每次迭代中，对于一组共同的粒子，同时在多个处理器上执行似然性评估，如算法1（第12行和第31行）所示。为了检验独立运行之间的变异性，展示了三个使用不同随机种子（SMC1–3）的独立估计；改变随机种子会影响算法1中的随机操作，包括粒子初始化和变异步骤（第4行、第29行和第34行）。同样，对于MCMC方法，展示了三个独立的估计（MCMC1–3），每个估计包含从不同初始值和随机种子生成的5000个样本。对于MCMC，采用了似然级并行性，其中对应于不同实验条件的似然贡献被并行评估，而马尔可夫链本身是顺序演化的。在这两种方法中，计算都使用了相同数量的计算核心（30个核心），以确保在计算努力方面的公平比较。MAP估计是通过使用差分进化算法最小化负对数后验来获得的，这是一种基于种群的全局随机优化方法。

3.2.2 使用实验数据的SMC方法估计
在案例2.2中，使用来自甲烷化过程的真实实验数据进行了参数估计。正如案例2.1中所展示的，由于参数之间的强相关性以及过程模型的高计算成本，这个估计问题具有挑战性，这可能会阻碍对后验分布的有效探索。在初步分析中，使用MCMC进行参数估计似乎只能获得有限的信息。然而，后续的研究表明，这种行为并不是MCMC方法的根本限制。当提供了适当的初始值时，Slice采样能够收敛到与使用SMC方法获得的后验分布基本一致的后验分布。这些结果表明，在适当的条件下，MCMC原则上可以解决这个问题。相应的结果在图S8和S9中报告。尽管如此，观察到MCMC的结果依赖于初始值的选择。对于MCMC3，在图6A中观察到后验分布明显变窄，这表明采样被限制在参数空间的一个局部区域。此外，尽管在图6C和D中，MCMC2的MAP估计在技术上被包含在后验分布中，但MAP估计附近的后验概率质量非常小，这意味着后验质量主要集中在远离MAP估计的区域。另一方面，在SMC方法中，真实值和MAP估计始终被包含在后验分布中，平均值的变化不像MCMC结果那样大。MAP估计始终位于后验密度高的区域这一事实表明，SMC方法能够以一种可靠的方式反映数据中包含的信息，而不依赖于初始值。从这些结果可以得出结论，SMC方法在独立运行中表现出更稳定的采样行为，而MCMC结果对随机种子和初始值的选择依赖性更强。在匹配的计算预算下，系统地、完全公平地比较SMC和MCMC之间的并行化策略仍然是未来工作的一个重要课题。相比之下，SMC框架中采用的并行化策略与MCMC中的链级并行化有根本不同。如算法1（第12行和第31行）所示，SMC在每次迭代中对一组共同的粒子并行化似然性评估。这种方法提高了计算效率，同时产生的结果与单核执行得到的结果相当。图5A展示了在单核条件下MH方法、Slice采样和SMC方法的平均计算时间比较。垂直轴上的平均时间代表使用不同随机种子进行的五次独立估计的平均计算时间，通过MH方法的平均计算时间进行了归一化。在这种归一化下，SMC方法的计算时间大约是MH方法的2.3倍，而Slice采样需要的时间大约是MH方法的4.3倍。

图5A展示了MH方法、Slice采样和SMC方法在单核条件下的平均计算时间比较。（A）计算时间比较，（B）通过并行实现减少SMC方法的计算时间。红色：SMC方法，蓝色：Slice采样，黑色：MH方法。图5B展示了并行化对SMC方法计算时间的影响，其中核心数量从1变化到60。计算时间大致与核心数量成反比减少，表明SMC方法具有高效的并行扩展能力。

3.2 甲烷化过程模型
关于第2.6.2节中描述的甲烷化过程模型的参数估计，我们比较了SMC方法和MCMC方法（Slice采样）。由于即使经过数十小时的估计，分布估计似乎也没有收敛到后验分布，因此没有考虑MH方法。

3.2.1 使用模拟数据的计算效率比较
案例2.1比较了每种方法的估计精度和计算时间，以定量评估计算效率的差异。SMC方法的超参数选择如下：，和。对于Slice采样，前2500个样本被丢弃作为烧录期。图6使用小提琴图展示了使用SMC和MCMC方法获得的采样结果。对于SMC方法，采用了粒子级并行性，即在每次迭代中，对于一组共同的粒子，同时在多个处理器上执行似然性评估，如算法1（第12行和第31行）所示。为了检验独立运行之间的变异性，展示了三个使用不同随机种子（SMC1–3）的独立估计；改变随机种子会影响算法1中的随机操作，包括粒子初始化和变异步骤（第4行、第29行和第34行）。同样，对于MCMC方法，展示了三个独立的估计（MCMC1–3），每个估计包含从不同初始值和随机种子生成的5000个样本。对于MCMC，采用了似然级并行性，其中对应于不同实验条件的似然贡献被并行评估，而马尔可夫链本身是顺序演化的。在这两种方法中，计算都使用了相同数量的计算核心（30个核心），以确保在计算努力方面的公平比较。MAP估计是通过使用差分进化算法最小化负对数后验来获得的，这是一种基于种群的全局随机优化方法。

图6展示了使用SMC方法和MCMC方法对案例2.1获得的采样结果进行比较。（A）、（B）、（C）和（D）。每个图中的水平线分别代表从上到下的第三四分位数、平均值和第一四分位数。SMC方法的结果用蓝色显示，MCMC方法的结果用红色显示。红色虚线代表真实值。粗绿色实线代表MAP估计，细绿色虚线代表伪数据的真实值。图6A–D显示了相似的整体趋势。对于MCMC1–3，真实参数值和MAP估计都包含在后验分布中。然而，不同估计运行之间观察到了显著的变异性，表明MCMC结果强烈依赖于初始值的选择。在图6A中观察到MCMC3的后验分布明显变窄，这表明采样被限制在参数空间的一个局部区域。此外，尽管在图6C和D中，MCMC2的MAP估计在技术上包含在后验分布中，但MAP估计附近的后验概率质量非常小，这意味着后验质量主要集中在远离MAP估计的区域。另一方面，在SMC方法中，真实值和MAP估计始终包含在后验分布中，平均值的变化不像MCMC结果那样大。MAP估计始终位于后验密度高的区域这一事实表明，SMC方法能够以一种可靠的方式反映数据中包含的信息，而不依赖于初始值。从这些结果可以得出结论，SMC方法在独立运行中实现了更稳定的采样行为，而MCMC结果对随机种子和初始值的选择依赖性更强。图S7提供了每种方法对实验误差σ的估计结果。在计算时间方面，观察到了明显的差异。SMC方法需要28小时进行估计，而MCMC方法大约需要150小时来获得5000个样本。因此，可以认为SMC方法在计算上是高效的，因为它能够在至少短一个数量级的时间内实现稳定采样。应该注意的是，由于实验噪声，后验分布与真实值之间的完全一致并不总是预期的。在比较数值推断方法的性能时，后验分布与真实值的精确收敛并不是主要关注点，因为这些值在实践中无法直接从噪声观测中获得。在这种情况下，将真实值包含在后验分布中主要是为了验证问题的表述是否在物理和统计上是一致的。相反，关键问题是推断方法是否能够可靠地识别出数据支持的高后验概率区域。从这个角度来看，MAP估计的位置特别具有信息性，因为它代表了被观测数据强烈支持的参数值。因此，在评估数值估计的质量和可靠性时，应该同时考虑后验分布的结构和MAP估计的位置。

3.2.2 使用实验数据的SMC方法估计
在案例2.2中，使用来自甲烷化过程的真实实验数据进行了参数估计。正如案例2.1中所展示的，由于参数之间的强相关性和过程模型的高计算成本，这个估计问题具有挑战性，这可能会阻碍对后验分布的有效探索。在初步分析中，使用MCMC进行参数估计似乎只能获得有限的信息。然而，后续的研究表明，这种行为并不是MCMC方法的根本限制。当提供了适当的初始值时，Slice采样能够收敛到与使用SMC方法获得的后验分布基本一致的后验分布。这些结果表明，在适当的条件下，MCMC原则上可以解决这个问题。相应的结果在图S8和S9中报告。尽管如此，观察到MCMC结果依赖于初始值的选择。对于某些初始化，采样轨迹倾向于局限于参数空间的有限区域，导致对后验分布的探索不足。这种行为在参数强相关和多模态结构的问题中很常见，这突出了在将MCMC方法应用于这类问题时初始化和调整的重要性。相比之下，SMC方法在所有测试场景中都表现出稳健和稳定的性能。如第3.2.1节中的图6所示，SMC结果对随机种子和初始粒子分布不敏感，始终产生解析度良好的后验估计。在案例2.2中，SMC从与MCMC相同的初始条件开始，能够在单次运行中完成参数估计。总计算时间大约为44小时，超参数设置如下：，和。这些选择遵循了第S5节中报告的敏感性分析，确保SMC设置不是为这个问题量身定制的。图7展示了参数估计的结果。如图7所示，除了某些参数外，后验分布的宽度比先验分布的宽度窄，这表明不确定性有所降低；也就是说，实验数据带来的信息增益得到了证实。此外，在不确定性降低的参数中，后验分布的均值与先验分布的均值有显著偏离，扩展到了先验分布中概率较小的区域。这表明所使用的实验数据包含了关于这些参数的强有力信息。除了图7中显示的边际分布外，还观察到几个参数对之间存在强相关性，反映了反应速率表达式的结构；详细的相关性模式在图S10中提供。

基于实验数据的估计。蓝色虚线代表先验分布，红色直方图代表后验分布。蓝色垂直线是先验分布的均值。(A)、(B)、(C)、(D)、(E)、(F)、(G)、(H)和(I)。

4 数据可用性和可重复性声明
所有支持本研究发现的数据都包含在附录S2中。用于生成图2-4（玩具问题示例）和图5-7（甲烷化案例）的数据也包含在这个档案中。与图3相关的核密度估计的额外结果也提供在内，这些结果出于简洁考虑在正文中被省略了。用于生成图S6的数据同样包含在附录S2中。所有数据集都是使用Python生成的。玩具问题示例也作为一个示例文件集在公共GitHub仓库中提供。对于甲烷化系统，如正文中提到的，模拟是使用assimulo库来求解微分-代数方程组的。微分方程和实验条件基于之前的研究。数据处理和图表生成的详细过程在正文中有所描述。

5 结论
在这项研究中，我们实现了一种结合了似然温度调整的SMC方法，并使用玩具问题和甲烷化过程模型进行了案例研究，以定量评估其性能，并将其与MCMC方法进行了比较。在使用玩具问题的案例研究中，我们使用SMC方法、MH方法和Slice采样从目标分布中抽取样本，并定量评估了分布的近似精度。结果表明，SMC方法在估计精度和计算时间上都更优。在使用甲烷化过程模型的案例研究中，使用模拟数据和实验数据对SMC方法和MCMC方法的计算效率进行了定量比较。在使用模拟数据的案例研究中，发现MCMC的估计结果随初始值的变化而大幅波动，并且需要较长的计算时间。另一方面，SMC估计可以在更短的时间内完成，从而证实了SMC方法的计算效率。在使用实验数据的案例研究中，观察到参数不确定性在先验分布和后验分布之间的降低和偏移，证实了实验数据提供的信息增益。此外，使用新操作条件下的实验数据验证了带有后验参数概率分布的模型。这项研究展示了SMC方法在参数相关性强和多模态估计问题中的实际优势，并表明该方法可以应用于真实的化学工程问题。同时，在似然评估计算成本较低且参数相关性较弱的情况下，MCMC方法可能更有效，表明SMC的相对优势取决于问题的特性。此外，还证明了SMC方法在计算效率上的优势。这一观察结果需要在未来的工作中通过其他化学过程进行验证。此外，系统地、公平地比较SMC和MCMC在相同计算预算下的并行化策略仍然是未来工作的一个重要课题。另外，由于实验数据较为有限，有必要使用具有不同操作条件和噪声水平的多个数据集来评估该方法。此外，还应与其他最先进的贝叶斯估计方法（如哈密顿蒙特卡洛和变分推断）进行比较。

作者贡献
Tatsuki Maruchi：软件开发；数据管理；调查；验证；形式分析；可视化；撰写——初稿。
Tomoyuki Yajima：软件开发；调查；验证；形式分析；监督。
Yoshiaki Kawajiri：概念化；撰写——审阅和编辑；方法论；调查；形式分析；监督；资金获取；资源管理；项目协调。

致谢
本研究部分得到了日本学术振兴会（JSPS）科学研究资助（项目编号JP25K01571）的支持。

数据可用性声明
支持本研究发现的数据可在以下链接公开获取：python-based-Sequential-Monte-Carlo-method-with-likelihood-tempering：https://github.com/maruchitatsuki/python-based-Sequential-Monte-Carlo-method-with-likelihood-tempering