计算昂贵的后验分布的高效 MCMC 方案

0. 摘要

我们考虑马尔可夫链蒙特卡洛 (MCMC) 计算方案，旨在最小化贝叶斯推理中后验分布的评估次数，当后验计算成本高昂时。我们的动机是计算昂贵的计算机模型的贝叶斯校准。先前在文献中提出的基于混合蒙特卡罗和目标分布的高斯过程近似的算法以三种方式扩展。首先，我们考虑将原始方法与回火方案相结合，以处理多峰后验分布。其次，我们考虑用高斯过程近似替换原始目标后验分布，这所需要的评估成本更低。第三，我们在回火方案的背景下考虑用高温链中的近似值替换真实目标分布，同时在最低温度链中保留真实目标。这在最低温度链中保留了正确的目标分布，同时避免了在涉及高温的移动中运行计算机模型的计算费用。我们的方法的应用被认为是一个降雨径流模型的校准，其中观察到参数后验的多模态。

1.介绍

环境科学和工程领域的研究人员广泛使用计算机模型来理解他们研究的系统并进行决策。这些计算机模型通常是确定性的，并在其实现中施加大量计算需求，通常需要数小时或数天才能运行一次。模型通常包含未知参数，需要根据不同模型参数值的计算机模型输出和观测数据的比较来估计这些参数：这意味着模型参数的校准在计算上要求很高，因为它需要重复运行电脑型号。

在统计框架中，通常计算机模型提供似然函数规范的一部分，而在贝叶斯校准方法中，将使用马尔可夫链蒙特卡洛 (MCMC) 方案来探索参数的后验分布。设计一个有效的 MCMC 方案以最小化计算机模型的评估次数是非常重要的，这也是激励我们工作的问题。然而，应该指出的是，本文中提出的想法更广泛地有用。我们的文章考虑了 Rasmussen (2003) 之前提出的方法的扩展。 Rasmussen 的想法是利用混合 Monte Carlo 提议（Duane et al. 1987）从目标后验分布进行模拟。混合蒙特卡洛需要重复评估目标分布对数的导数来进行建议移动，并且拉斯穆森考虑在生成建议时用高斯过程近似替换真实目标。在蒙特卡洛方案的接受/拒绝步骤中，使用真正的目标分布，保持正确的遍历分布。

我们的文章将 Rasmussen 的工作扩展到参数分布是多模态的情况，这是他不打算解决的情况。我们通过使用高斯过程近似与平行回火结合混合蒙特卡罗的想法来做到这一点（Geyer 1991；Hukushima 和 Nemoto 1996）。 Chen, Qin, and Liu (2000) 也考虑了混合蒙特卡罗与平行回火相结合，但没有使用高斯过程近似。在平行回火中，许多链条平行运行，由一个在物理类比中通常称为温度的参数索引。不同链的目标分布是原始目标分布的缓和或过度分散版本，最低温度链对应于感兴趣的原始目标分布。每条链都可以单独更新，但也可以考虑“交换”动作，为不同的链交换状态。这意味着可以在过度分散的高温链中迈出大步，然后这些状态可以向下交换到最低温度链。这个想法允许我们在原始目标分布的孤立模式之间移动。通常，计算机模型的校准问题会导致复杂的多峰后验分布——稍后我们将研究这种情况下的降雨径流模型的应用。我们在本文中的第一个贡献是 Rasmussen 的想法的基本扩展（扩展到参数分布是多模态的情况），包括平行回火。一个关键的观察结果是，可以在所有温度下重复使用对回火目标密度之一的评估，以构建对所有温度下的提议有用的高斯过程近似。我们的第二个贡献是检查一个想法（类似于 Bliznyuk 等人 2008 的想法），用高斯过程近似代替真实目标分布以进行采样，而不仅仅是生成建议移动。我们的第三个贡献是在回火方案中应用这种替换真实目标的想法。然而，我们建议在最低温度链中保留真实目标（这是感兴趣的后验分布）。这在最低温度链中保持正确的目标分布，但避免了在涉及最高温度链的移动中运行计算机模型的费用。这个想法非常重要，因为这意味着如果计算复杂度由计算机模型的评估次数支配，那么这个版本的回火计算成本只有 Rasmussen 原始方案仅涉及单个链的计算成本的两倍左右。

许多其他作者已经考虑过计算机模型的贝叶斯校准。与我们最密切相关的工作是 Bliznyuk 等人的工作（2008 年）。他们考虑对感兴趣的目标分布进行径向基函数逼近，并在参数空间中明确定义设计区域以构建逼近。他们认为，如果运行计算机模型的计算预算是固定的，那么所有预算都应该分配给近似值的构建，然后采样应该基于对目标分布的近似值进行。如果计算预算真的是固定的，我们会同意，但在实践中，有时可以灵活地根据需要分配更多的计算资源，我们相信 Rasmussen (2003) 使用真实目标的方法是有帮助的，因为采样阶段可以帮助诊断质量差的近似值。布利兹纽克等人(2008) 提出了一些额外的新想法，例如区分计算机模型参数和其他在更改时不需要重新运行计算机模型的参数，以及通过剖析降低维度的各种策略。降维很重要，因为对目标分布的任何非参数近似只能在相当低的维度（比如小于 10 维）中很好地工作。这些降维策略也可以与我们文章中描述的方法一起应用。 Bliznyuk 等人的方法。 (2008) 不打算与多模式目标分布一起使用。

在其他关于计算机模型贝叶斯校准的工作中，Kennedy 和 O’Hagan (2001) 制定的框架特别有影响力。他们建立在使用高斯过程对计算机代码输出建模的想法（Sacks et al. 1989）。 Kennedy 和 O’Hagan (2001) 的公式允许各种不确定性来源，包括但不限于模型输入中的不确定性、测量误差、代码不确定性和物理系统的不充分表示。一个关键的想法是使用模拟器处理代码不确定性（这是由于只能对计算机模型进行有限数量的评估而导致的不确定性），这是对模拟模型的统计近似，可以以更少的计算成本进行评估.在这里（与 Rasmussen 2003 和 Bliznyuk 等人 2008 类似）我们实际上使用仿真器来描述后验分布本身，而不是模拟输出。根据可用数据和所研究物理系统的性质，在 Kennedy 和 O’Hagan (2001) 的框架中考虑的所有不确定性来源的单独建模可能是不可能或没有必要的。在我们后面的真实示例中，我们使用类似于 Bliznyuk 等人的非线性回归公式。 (2008) 具有相关误差，其中误差代表测量误差、模型不足和观测输入测量误差的组合。肯尼迪和奥哈根（2001）的基本框架已被广泛使用并以各种方式扩展。参见 Higdon 等人。 （2008 年）回顾了当前的技术状态并将基本框架扩展到涉及高维输出的情况。 Higdon, Lee, and Holloman (2003) 描述了另一种与我们密切相关的贝叶斯校准方法。这里考虑了在计算机模型的贝叶斯校准中运行并行 MCMC 链的想法，其中并行链源自计算机模型在不同分辨率级别的运行以及相应的不同计算费用级别。
在下一节中，我们将回顾混合蒙特卡罗和平行回火，这是我们建议的方法的基础。第 3 节总结了 Rasmussen (2003) 的方法。第 4 节描述了我们的扩展，第 5 节考虑了一些模拟示例，第 6 节讨论了在后验分布的多模态明显的降雨径流模型校准中的应用。第 7 节结束。实现本文中讨论的算法的 R 包可在 http://cran.r-project.org/web/packages/MCMChybridGP 获得。

2.马尔科夫链蒙特卡罗

马尔可夫链蒙特卡罗是一种从感兴趣的分布中生成样本的方法。 在贝叶斯推理应用中，这是一组未知数 $\theta$ 上的后验分布 $p(\theta |y)$，例如校准练习中的模型参数。 MCMC 在参数空间 { θ ( n ) : n ≥ 0 } \{θ(n) : n ≥ 0\} {θ(n):n≥0} 上构造马尔可夫链，使得链的平稳分布为 $p(\theta |y)$。 然后通过采用任意起始状态 ${\theta }^{(0)}$ 并长时间运行链，并在丢弃样本不是典型的平稳分布的“老化”期后，我们从 $p(\theta |y)$获得依赖样本 . 这里假设您基本熟悉 MCMC。 例如，参见 Liu (2001) 的介绍。

2.1 混合蒙特卡罗

几乎所有在实践中使用的 MCMC 算法都是 Metropolis-Hastings (M-H) 算法的变体，混合蒙特卡洛涉及在增强状态空间上考虑人工动力系统，作为在 M-H 算法中生成提议移动的一种方式 . 考虑一个分布 $\pi (x),x=(x_1,...,x_n)$。 下面我们假设 $\pi (x)$ 中没有未知的归一化常数，但如果不是这种情况，那么下面的讨论几乎没有变化。 将 $\pi (x)$ 写为 $\pi (x)=exp(-Z(x))$ 其中 $Z(x)=-log\pi (x)$。 在动力系统类比中，$Z(x)$ 称为势能。 考虑一组独立于 x 的辅助变量 $p=(p1,...,pn)$，分布

和联合分布 $\pi (x,p)=\pi (x)\pi (p)$。 变量 $p$ 称为动量，$W(p)={\sum }_{i=1}^n{p}_i^2/2$ 是动能。 那么 $\pi (x,p)=(2\pi {)}^{-n/2}exp(-H(x,p))$ 其中 $H(x,p)=Z(x)+W(p)$ 并且$H(x,p)$ 是 称为哈密顿量。 考虑根据哈密顿方程演化的状态 $(x,p)$ 的人工动力系统

其中 $t$ 是人为的时间参数。 根据 (1) 的状态演化使状态的分布 $\pi (x,p)$ 不变。 该性质与哈密顿动力学的能量保持、体积保持和时间可逆性有关。 参见 Liu（2001，第 9 章），从统计角度介绍混合蒙特卡洛，包括对牛顿力学的简要介绍。 另见尼尔 (2010)。 如果我们有一个来自 $\pi (x,p)$的样本，那么 x 的边际分布是 $\pi (x)$，所以我们可以丢弃 p 的模拟值以从 $\pi (x)$ 中获得一个样本。

通常我们不能解析解方程（1），我们必须用数值近似。 可以使用所谓的跳跃式离散化来做到这一点并保留时间可逆性。 假设我们知道在时间 t 的状态$(x,p)$，并且我们想要根据 (1) 来近似在时间 $t+\delta$ 演化的状态。 在越级离散化中，我们允许 $x_i(t+\delta )$ 和 $p_i(t+\delta )$ 近似为

随着 $\delta$ 趋于零，根据 (2) 演化的 $x_i(t)$ 和 $p_i(t)$ 将匹配根据 (1) 演化的轨迹。然而，对于有限 $\delta$，情况并非如此，根据离散动力学 (2) 演化的轨迹不会使 $\pi (x,p)$ 保持不变。 Duane 等人的混合蒙特卡罗的想法。 (1987) 是使用 Metropolis-Hastings 接受/拒绝步骤来纠正这一点，以保持 $\pi (x,p)$ 作为不变分布：在每次迭代中，当前值 $(x^C,p^C)$ 更新如下: 1. 通过对 $L$ 个大小为 $s\delta$ 的步长使用跳跃式离散化生成建议状态 $(x^P,p^P)$，其中 $\delta >0$ 是步长长度，$s$ 选择为 1 或 -1，每个概率为 0.5。 2. 以概率 m i n { 1 , α } min\{1,α\} min{1,α}$ 接受 $(x^P,p^P)$ 作为 $(x,p)$ 的新值，其中 α = e x p { − H ( x P , p P ) + H ( x C , p C ) } α = exp\{−H(x^P, p^P) + H(x^C, p^C)\} α=exp{−H(xP,pP)+H(xC,pC)}。否则 $(x,p)$ 的值保持在 $(x^C,p^C)$。在这里，Lapfrog 离散化的大小为 $\delta$ 的步长用于获得 Metropolis-Hastings 算法的建议。 $L$ 和 $\delta$ 需要由用户选择。稍后我们将固定 L 并通过反复试验选择 δ 以达到一定的提案接受率。在计算中，如果我们考虑越级算法的多个步骤，则 (2) 中指示的 p 的两个半步可以合并为一个。

2.2 并行回火

在传统的 MCMC 中，构建了一个马尔可夫链，其中目标分布是后验分布。 写为：

其中 $p(\theta )$ 是 $\theta$ 的先验分布，$p(y|\theta )$ 是可能性。 传统的 MCMC 算法如果由于难以构建好的提议而具有多个孤立模式，则难以探索目标分布。 平行回火（Geyer 1991；Hukushima 和 Nemoto 1996）考虑了一系列分布

$j=1,...,N$，其中 $1=T_1<\cdot \cdot \cdot <T_N$ 是一组所谓的温度参数。 在高温下，${\pi }^j(\theta )$ 接近均匀，$p(\theta |y)$ 中的任何孤立模式都被展平，而 ${\pi }^1(\theta )$ 在 $T_1=1$ 时为 $p(\theta |y)$。 如果$\Theta$ 表示参数空间， 考虑乘积空间 ${\Theta }^N$ 上的分布 $\pi (\cdot )$，其中让 ${\theta }^S=({\theta }^1,...,{\theta }^N)\in {\Theta }^N$ 且每个 ${\theta }^j\in \Theta$

我们可以通过 MCMC 从$\pi ({\theta }^S)$ 生成样本。 如果获得来自 $\pi ({\theta }^S)$ 的样本，则保留模拟的 ${\theta }^1$ 值并丢弃 ${\theta }^S$ 的其余部分会得到来自 $p(\theta |y)$ 的样本。 在并行回火算法的最常见形式中，在随机选择其中一个 ${\theta }^j$ 以通过传统的 M-H 步骤更新与 ${\theta }^S$ 的其余分量保持固定之间进行交替，以及提议交换状态以对应于相邻的分布的移动 温度参数（根据通常的 M-H 规则接受/拒绝交换）。 高温下的大变化可以通过交换移动传递到低温，从而允许在隔离模式之间移动。 通常选择间隔的温度，以便 $T_{j+1}/T_j=c>0$ ，对于某个值 c$j=1，.....,N-1$ 。

3. 拉斯穆森法

混合蒙特卡罗方法需要计算目标密度对数的导数。 这在许多应用中可能是可行的，或者至少在计算上要求很高，Rasmussen (2003) 提出了使用对数目标密度的高斯过程近似并使用近似的导数在混合蒙特卡罗中构建提议的想法。 我们在这里不介绍高斯过程，但请读者参考 Rasmussen 和Williams (2006) 的背景知识。 在 Rasmussen 的方法中，使用高斯过程来描述关于对数目标密度 $f(\theta )$ 的不确定性。 对数目标密度已在某些点进行评估，给出 $y_i=f({\theta }_i),i=1,..,n.$ 我们使用 $f(\theta )$ 的先验，它是一个高斯过程，平均 $E(f(\theta ))=0$ 和一个平稳的协方差函数 $v(h)$，使得

其中 ${\lambda }^2$ 是方差参数，$r(\cdot ,\eta )$ 是取决于参数 η 的相关函数。 我们稍后在示例中选择 $r$ 作为Matérn 类的一个特例。

在观察到 $y=(f({\theta }_1),...,f({\theta }_n))$ 之后，如果已知 ${\lambda }^2$ 和 $\eta$，则函数 $f$ 的后验分布也是一个高斯过程。我们将 $f$ 的后验均值写为 $\mu (\theta )=E(f(\theta )|y)$，后验标准差函数写为 $\sigma (\theta )=\sqrt{Var(f(\theta )|y)}$。当然，可以将先验放在未知参数 ${\lambda }^2$ 和 $\eta$ 上，并考虑更复杂的分层模型结构，但我们在这里不这样做（参见，例如，Rasmussen 和 Williams 2006，进一步讨论）。我们需要某种方法来获得 ${\lambda }^2$ 和 $\eta$ 的点估计。在 Rasmussen (2003) 之后，这是通过边际最大似然来完成的——边际对数似然及其导数很容易计算，并且可以使用标准优化技术。人们还可以考虑观察到 $f$ 的导数的情况，因为 $f$ 及其导数（假设它们存在）之间的互协方差函数很容易为高斯过程计算。在使用导数信息逼近 MCMC 目标分布的应用中，尤其是在高维中，这种扩展可能是有益的，Rasmussen (2003) 考虑了这种扩展。然而，对于我们以后的降雨径流应用，目标密度导数的评估是不可行的，我们不进一步考虑这一点。 Rasmussen 的算法由两个步骤组成：探索阶段和采样阶段。在探索阶段，我们执行以下步骤。

算法1：

1. 生成一组初始随机点 ${\theta }_1,..,{\theta }_i$。这些可能是随机生成的（例如从先前的），也可能是确定性地选择的。
2. 评估初始点处目标密度的对数（用 $\pi (\theta )$ 表示目标密度，可能达到未知的乘法常数），$y_j=\log \pi ({\theta }_j),j=1,..,i$
3. 对点 $(y_j^{\ast },{\theta }^j)$ 拟合高斯过程，其中 $y_j^{\ast }=y_j-\overline{y}$ 和 $\overline{y}=n^{-1}{\sum }_{k=1}^i{y}_k$。
4. 将 ${\theta }^C$ 初始化为步骤 1 中的状态之一，随机均匀选择。
5. 对于 $j=1，...,E$（其中 $E$ 表示探索阶段样本大小），使用从 ${\theta }^C$ 开始的 $L$ 步的跳蛙离散化（在 Rasmussen 之后，我们在示例中使用 L = 1000），步长为 $\delta$ ，以获得建议状态 ${\theta }^P$。建议模拟中使用的目标分布不是 $\pi (\theta )$，而是与 $\exp (\mu (\theta )+s_p\sigma (\theta )/\lambda )$ 成正比，其中 $\mu (\theta )=E(f(\theta )|y)$, $\sigma (\theta )=\sqrt{Var(f(\theta )|y)}$，$f$ 表示我们对 $\log \pi (\theta )$ 的高斯过程近似，$y$ 表示到目前为止对目标分布的评估，$\lambda$是高斯过程协方差中 $\lambda$ 的当前估计值函数，$s_p$ 是一个调整参数（默认值 1），可用于控制提案移动的分散程度。在每个越级步之后监控状态，如果 $\sigma (\theta )/\lambda >0.9$，我们在全部步数之前停止，其中 $\theta$ 表示当前状态。在最终状态 ${\theta }^P$ 我们评估 $\pi ({\theta }^P)$，更新高斯过程近似，并更新高斯过程近似中的协方差函数参数。然后我们以概率 m i n { 1 , α } min\{1,α\} min{1,α} 将 ${\theta }^C$ 更新为 ${\theta }^P$，其中
   
   与 Rasmussen 相比，上述算法存在一些细微差别。这些与我们使用 $s_p\sigma (\theta )/\lambda$ 而不是 $\sigma (\theta )$ 的目标分布以及在越级步骤中提前停止的相应标准有关。 Rasmussen 在探索阶段引入了修改后的目标分布，以鼓励在目标分布值的不确定性很大的区域中探索空间。允许一个可以修改这种效果强度的额外调整参数肯定会使我们的方法更适应手头的问题。然而，这种适应性可能是喜忧参半，因为需要设置此调整参数。拉斯穆森最初的提议是一个明智的默认。进一步说明，在步骤 1 确定性地选择初始点集通常也很有帮助。事实上，在这种情况下，我可能会被选为探索阶段考虑的点总数的很大一部分。在后验分布具有单一模式的情况下，Bliznyuk 等人。 （2008 年）当目标是使用灵活的回归方法近似后验分布时，考虑对评估先验时间似然性的点进行确定性设计。另一种方法是选择一组初始点来广泛覆盖空间（可能使用拉丁超立方体设计，Morris 和 Mitchell 1995），然后通过 Rasmussen 探索阶段的随机搜索来细化这组点。一旦探索步骤完成，下一阶段是采样阶段，其进行如下。

算法2:

1. 将 $\theta ={\theta }^C={\theta }^{(1)}$ 初始化，利用最大化 $\log \pi (\theta )$，在探索阶段找到的 $\theta$ 值。
2. 对于 $j=2，...,s$
   • (a) 根据 $\mu (\theta )=E(f(\theta )|y)$ 对 $\log \pi (\theta )$ 的逼近，通过从 ${\theta }^C$ 开始的跳跃式离散化，以步长 $\delta$ 的 $L$ 步生成建议值 ${\theta }^P$ .
   • (b) 以概率 m i n { 1 , α } min\{1,α\} min{1,α} 接受 ${\theta }^{(j)}={\theta }^C={\theta }^P$，其中 $\alpha =\pi ({\theta }^P)/\pi ({\theta }^C)$ ，否则 ${\theta }^{(j)}={\theta }^C={\theta }^{(j-1)}$。

该算法的最终结果是一组值 ${\theta }^{(1)},..,{\theta }^{(s)}$ 是来自 $\pi (\theta )$ 的近似样本。

上述方法的一个限制是在对数目标分布的高斯过程近似中使用平稳协方差函数。 虽然结合问题知识的更复杂的协方差模型可能会带来好处，但 Rasmussen (2003) 观察到，在高密度区域中，静止模型可能是合理的。我们简要讨论第 6 节示例中出现的另一个问题。 必须从参数值有界的后验分布进行模拟，并且通常在混合蒙特卡罗中假设无界空间。 最近 Neal (2010) 概述了一种在这种情况下实现混合蒙特卡罗的技术，其中约束被视为反映障碍，并且可以直接将我们的算法与 Neal 的方法结合使用。

4. 扩展

我们对 Rasmussen 算法的扩展分为三种，在下面的三个小节中进行解释

4.1 并行回火和混合蒙特卡罗

首先，我们考虑将混合蒙特卡罗和高斯过程近似与并行回火相结合，以给出一种适用于可能的多模态目标分布的方法。我们按如下方式实现。在算法 1 的第 2 步中，我们按照原始算法中的方法继续构建对 $p^1(\theta )$ 的高斯过程逼近，其中 $p^k(\theta ),k=1,..,N$, 是第 2.2 节中描述的回火分布。请注意，$\log p^k(\theta )=1/T_k\log p^1(\theta )$，因此对 $\log p^1(\theta )$ 的高斯过程近似立即给出对 $\log p^k(\theta )$ 的近似。这是我们在并行回火算法中用于在不同温度下更新链的越级步骤的近似值。在算法 1 的第 4 步中，我们需要每个温度下的状态初始值。我们以与 { p ( θ j ) p ( y ∣ θ j ) } 1 / T k \{p(θ_j)p(y|θ_j)\}^{1/T_k} {p(θj​)p(y∣θj​)}1/Tk​ 成比例的概率从步骤 1 的状态中随机采样，以获得温度 $k$ 的初始值。在算法 1 的第 5 步中，我们现在必须使用 L 步跨越式离散化进行 N 次单独更新，每个温度更新一次。在每一步之后，我们根据对所有温度下目标分布的评估，将高斯过程近似更新为 $\log p^1(\theta )$。请注意，如果我们在某个 $\theta$ 值处评估 $p^k(\theta ),k>1$，我们已经在相同 $\theta$ 处评估了 $p^1(\theta )$，我们可以相应地更新我们的高斯过程近似为 $\log p^1(\theta )$。我们还在步骤 5 中考虑了不同温度下的状态交换步骤。可能会随机选择一对相邻的温度进行交换，或者可能在每个步骤中考虑所有相邻温度对的交换。在算法 2 的第 2 步中对每个温度和交换移动进行了类似的更新。

4.2 近似抽样

在算法 2 的第 2 步中，我们可以用高斯过程近似替换真实的目标分布[即，我们可以在接受概率中使用 $exp(\mu (\theta ))$ 代替 $\pi (\theta )$]。 Bliznyuk 等人也考虑了类似的想法。 （2008）他们考虑了目标分布对数的径向基函数近似，尽管他们的方法不适合多模态目标分布。当 $p(y|\theta )$ 的评估代价高昂时，用近似值替换目标分布的想法可以节省计算量，代价是仅从目标近似地生成样本。布利兹纽克等人(2008) 指出，近似需要限制在对近似质量有信心的空间区域。在这里，我们不对支持施加硬约束，而是修改了我们使用高斯过程近似来构造适当目标分布的方式，以便我们不鼓励在目标分布不确定性很大的部分空间中进行采样。我们考虑未归一化的目标密度 $exp(\nu (\theta ))$，其中

我们使用符号 z ＋ = m a x { 0 , z } z_＋ = max\{0, z\} z＋​=max{0,z} 表示实数$z$ 的正数部分，我们使用 $2$ 作为 $r$ 的默认值。

请注意，如果 $\sigma (\theta )/\lambda <0.9$，我们的新目标分布会减少到 $exp(\mu (\theta ))$，但是会增加一个惩罚，该惩罚会随着 $\sigma (\theta )/\lambda$ 接近 1 而增加，以阻止在空间部分进行采样，其中 高斯过程近似的方差很大。 显然，0.9 的阈值可以在上述结构中变化。

4.3 仅在最低温度下更换真实目标

在结合我们的前两个回火修改和替换真实目标时，我们进一步建议保留最低温度链的真实目标，同时替换高温链中的近似值。这意味着只有在涉及最低温度的更新（单独更新到最低温度状态或涉及最低温度状态的交换）时，我们才需要运行计算机模型。如果计算机模型运行的计算时间占总计算时间的主导，那么这个版本的回火计算复杂度大约是运行单链的两倍。此外，由于我们保留了最低温度的真实目标，因此我们准确地为这个最低温度链保持了感兴趣的平稳分布。但是请注意，这种修改对于发现在探索阶段根本没有访问过的模式并不是特别有用。这样做的原因是我们依靠交换步骤在不同模式之间移动，并且仅当最低温度的状态对应于替换目标中的一种模式（即，与在探索阶段发现的模式）。

为了实现我们的修改，在算法 2 的步骤 2(b) 中，我们在采样阶段使用 M-H 接受概率

对于最低温度链的单独更新，对于温度 $T_k，k>1$ 的单独更新，

其中 $\nu (\theta )$ 在第 4.2 节中定义。 对于平行回火交换步骤 α 变为

建议交换最低两个温度链的当前值 $θ_1 和 ${\theta }_2$。 对于涉及温度 $T_j$ 和 $T_{j+1}，j>1$的交换，α 变为

5. 模拟示例

我们考虑两个有启发性的模拟示例。 第一个考虑七个维度的多峰目标分布，第二个考虑一个单峰目标，但在目标密度具有扭曲的非椭圆轮廓的意义上，依赖结构非常复杂。 类似于最后一种类型的目标密度经常出现在非线性回归问题中。 我们在下面报告仅运行单个链的结果，以显示单个运行中的混合，但性能类似于跨不同复制报告的结果。 为了便于与其他方法进行比较，可以使用实现我们方法的 R 包（http://cran.r-project.org/web/packages/MCMChybridGP）。

5.1 七个维度的多模式示例

为了说明我们在多模态情况下的方法，我们考虑从七维正态分布的两分量混合进行模拟的问题。 特别是，我们让

其中 $p_1=p_2=0.5,{\mu }_1=(1,1,0,0,0,0,0),{\mu }_2=-{\mu }_1,{\Sigma }_1={\Sigma }_2=0.35I_7$ 其中 $I_7$ 表示 7 × 7 单位矩阵.

我们考虑了四种不同的采样算法：
Rasmussen 的原始方法（方法 M1） Rasmussen 算法的并行回火版本，真实目标用于所有温度下的接受/拒绝步骤（方法 M2） 具有高斯的并行回火方法在所有温度下的接受/拒绝步骤中使用的过程近似（方法 M3）和使用高斯过程近似的平行回火方法用于除最低温度链之外的所有目标（方法 M4）。在比较方案时，我们在方法 M2 的探索阶段使用了 200 次迭代，然后使用为该方法学习的高斯过程近似作为所有方法的采样阶段的近似。这样做的原因是，我们想在所有方法的探索阶段获得的信息质量相同的情况下，比较这些方案。在回火方案中，使用了五个温度，$T_{j+1}/T_j=1.5，j=1,..,4$. 在跳跃步骤中，我们使用 $L=1000$ 和 $\delta =0.001$。在算法的初始化中，我们使用了从最大拉丁超立方体设计中生成的 20 个点。对于每个变量，生成的初始点位于区间 [-2, 2] 内。图 1 是 5000 次采样迭代和方法 M1、M2、M3 和 M4 的 $x_1$ 和 $x_2$ 与迭代的轨迹图。自相关函数如图 2 所示。

图 1. 方法 M1（左上）、M2（右上）、M3（左下）和 M4（右下）在采样阶段 x1 和 x2 的迭代图。

图 2. 方法 M1（左上）、M2（右上）、M3（左下）和 M4（右下）在采样阶段迭代 x1 和 x2 的自相关函数图。 该图的在线版本是彩色的。

图 1 中的主要观察结果是方法 M1 在两种模式之间不能很好地混合，这是 Rasmussen (2003) 的原始算法没有做的事情，也不是为之设计的。这样做的一个结果是，在我们 5000 次迭代的计算“预算”内，方法 M1 在两种不同模式中花费的时间可能无法准确反映两种模式的概率质量。图 2 的自相关函数，特别是它们对于不同方法衰减到零的速度有多快，可以提供有关不同方法的后验不相关样本有多接近以及估计感兴趣的后验期望的精度的信息。显然，方法 M2-M4 远远优于方法 M1。我们的扩展方法很好地表示了目标分布，但是在所有温度下都使用高斯过程近似的方法 M3 中显示了一些准确性损失——这从查看边缘分布的非参数密度估计中可以明显看出（结果未显示）。方法 M4 在混合方面几乎与 M2 一样好，但是就真实目标分布的评估次数而言，方法 M4 的计算开销当然要少得多。

5.2 一个非椭圆轮廓的例子

接下来，我们考虑具有复杂依赖性的二维密度，其中密度的轮廓是非椭圆的。我们的示例取自 Haario、Saksman 和 Tamminen (2001)，他们描述了一种将正常密度转换为香蕉形轮廓的方法。构造通过取一个多元正态密度 $h(x)$ 来进行，比如说零均值和协方差矩阵 $C=diag(100,1,...,1)$，然后对于参数 $b>0$，我们可以定义$h_b(x)=h◦\varphi b$ 其中 ${\varphi }_b(x)=(x_1,x_2+b{x}_1^2-100b,x_3,...,x_p{)}^T$。这里 $h_b(x)$ 是一个密度函数，参数 $b$ 控制轮廓的扭曲程度，$b=0$ 对应于正态性。这里我们只考虑 $p=2$ 和 $b=0.03$ 的情况。我们使用 5000 次采样迭代和方法 M1、M2、M3 和 M4 以及与前两个示例相同的算法参数设置，在探索阶段进行 1000 次迭代。

迭代的散点图如图 3 所示，自相关函数如图 4 所示。

• 图 3. 方法 M1（左上）、M2（右上）、M3（左下）和 M4（右下）在采样阶段 x1 和 x2 的迭代图。

从图 3 中我们可以看出，所有方法都能够探索目标分布，尽管在方法 M3 的情况下，目标被替换为尾部没有被充分探索。与 Rasmussen 的原始方案（方法 M1）相比，方法 M2 和 M4 都混合得很好，但和以前一样，方法 M4 比 M2 涉及更少的计算开销。

• 图 4. 方法 M1（左上）、M2（右上）、M3（左下）和 M4（右下）在采样阶段迭代 x1 和 x2 的自相关函数图。 该图的在线版本是彩色的。

图 4 显示，与方法 M1 相比，方案 M2-M4 的 MCMC 迭代自相关衰减到零的速度要快得多。这意味着对于方法 M2-M4 的给定数量的 MCMC 迭代，对感兴趣的后验期望的估计将更加精确。

6. 应用：降雨-径流建模

我们描述了我们的方法在降雨径流模型校准中的应用。水文学家使用降雨径流模型来描述集水区对降雨事件的响应。在这里，我们关注 SOBEK-Urban 径流模型在新加坡单个子汇水面积中的应用。参见 Verwey 等人 (2008) 进一步讨论该模型及其在新加坡背景下的应用。我们关注的是模拟子汇水面积中单个主河道中的流量以响应降雨事件。 SOBEK-Urban 使用“沙井”和“管道”网络对降雨径流过程进行了概念化。沙井储存总径流（或作为径流的收集点），管道将储存的径流输送到主渠道。图 5(a) 显示了一个带有排水网络的子汇水面积，它被图 5(b) 中概念化的人孔和管道系统所取代。子汇水面积表面到人孔的径流由公式 $q=ch$ 描述，其中 $q$ 是从汇水面积进入人孔的流入量（mm/min），c 是线性径流因子（$mi{n}^{-1}$），$h$ 是集水区的蓄水量（毫米）。

径流因子 $c$ 是流域特征（长度、粗糙度和坡度）的函数，基于这些 SOBEK-Urban 有 12 种具有不同 c 值的区域类型。在这个例子中，只使用了两种面积类型，给出了两个径流因子，它们的值需要校准。线性地表径流首先储存在人孔中，然后通过管道进入主渠道。定义沙井的参数是街道层的存储区域、底层的存储区域和沙井的底层（或沙井内底水平）。表 1 总结了模型中的参数及其范围。校准中这些参数的先验分布被认为在给定范围内是均匀的。图 6 显示了校准数据。

数据显示了 2005 年 12 月 17 日至 20 日期间，新加坡滨海湾每 10 分钟的流量。这一时期对应于一次重大的洪水事件。请注意，流量可以是负数——这是因为在潮汐河流中，流量既可以流入大海，也可以流出大海。我们的数据模型采用以下形式

其中 $Q_t$ 是时间 $t$ 的流量，$f(z_t;\theta )$ 是时间 t 的降雨径流模型输出，$z_t$ 表示降雨输入，$\theta$ 表示模型参数。 两侧的变换是为了更接近误差 $t$ 的正态性和同方差性，假设在 $Cov(t,t)=exp(-a|t-t^{\prime }|)>0$。我们在 ${\sigma }^2$, $IG(0.01,0.01)$ 上使用共轭逆伽马先验，它允许 ${\sigma }^2$ 从后验分析中积分。 这里有大量关于 ${\sigma }^2$ 的可能性信息，因此关于这个参数的推断不会对分析方便的模糊共轭先验规范敏感。 我们还修复了一个点估计，给出了大约 0.1 的滞后一自相关。 这是基于使用先前相关应用的校准对模型拟合的残差进行的分析。 我们专注于探索七个计算机模型参数的后验分布。

这个例子是一个非常具有挑战性的例子。即使找到模型参数的点估计作为似然面的后验模式或最大值也不是直截了当的，因为似然面是多模态的。多模态优化问题通常使用类似于我们在这里描述的随机搜索方法来解决，但适用于优化而不是后验模拟。我们探索后验的目标不仅仅是发现似然或边际后验的主要模式，而是根据它们的后验概率以正确的比例对各种模式进行采样。因此，即使在算法的探索阶段发现了模式，如果我们想要探索复杂的后验分布和似然面，后续的采样也是必要的。对于算法的初始化，我们通过最大拉丁超立方体设计生成 100 个点。我们使用 L = 1000 和 δ = 0.00001。使用了 400 次迭代的探索阶段，然后运行了 25,000 次采样迭代，前 5000 次被丢弃为“老化”。方法 M1 和 M4 的所有变量的轨迹图如图 7 所示，

这表明方法 M4 的交换移动能够在几种不同的模式之间快速切换，而方法 M1 不能做到这一点（例如，参见轨迹图对于${\theta }_6$）。检查几个变量的联合变化显示了复杂的多峰结构。我们仅展示了其中两个参数的一些结果：管道宽度的比例因子 (${\theta }_4$) 和渗透区域的径流因子 (${\theta }_7$)。图 8 显示了所有四种不同方案的 $({\theta }_4,{\theta }_7)$ 迭代的散点图。

Sobek 模型的每次运行大约需要 15 秒，对于我们计算成本最高的方法 (M2)，25,000 次采样迭代的总计算时间约为 14 天。在我们 25,000 次迭代的计算预算内，是否有任何方法充分探索了这里的后验分布，这是值得怀疑的。但是，很明显，仅涉及运行单个链的方法 M1 是无效的。令人怀疑的是，方法 M3 的“替换目标”是否代表这个复杂的 7 维后验分布的后验井，探索阶段为 400 次迭代。方法 M2 和 M4 探索了几种不同的模式。

图 9 显示了 $({\theta }_4,{\theta }_7)$ 的所有四种不同方案的迭代轨迹图，显然方法 M2 比方法 M4 混合更快，但方法 M2 涉及评估计算机模型在所有温度和运行的计算成本要高得多。图 10 还显示了不同方案的已访问状态的对数似然。

对于方法 M1、M2 和 M4，图表在 y 轴上使用相同的比例，但对于方法 M3，比例不同，因为对数似然值的范围非常不同（方法 M3 访问后验概率非常低的区域，因为用近似值替换真实目标）。请注意，方法 M3 不会计算真实目标的实际对数似然，因此这些值是在采样完成后计算的。方法 M2 和 M4 是比方法 M1 具有更高对数似然的访问模式。

模型校准参数的估计对于任何流域都非常重要。对于集中时间非常短的流域（即流量传播时间；例如新加坡的 30 分钟），更重要的是模型能够以足够高的准确度预测峰值流量/流量及其到达时间，将提供有用的信息，例如洪水预警措施。我们在这个例子中描述的用于校准和探索复杂多模态后验的方法在这种情况下非常重要。

7. 结论和未来工作

已经提出了 Rasmussen (2003) 方法的扩展，该方法结合了混合蒙特卡罗、目标密度的高斯过程近似和平行回火。该扩展在要采样的目标密度是多峰的问题中很有用。还有一些剩余的挑战为未来的工作提供了机会。特别是，看看该方法的基本思想是否可以在高维问题中更有用将会很有趣。与其使用一般的高斯过程方法来近似对数目标密度，不如在参数数量较大时，对目标密度进行更结构化的近似可以在混合蒙特卡罗机制中提供有用的建议。 Rasmussen (2003) 也提出了在高斯过程近似中包含更多问题特定知识的静态协方差结构的替代方案也可能是有益的。