面向分析的随机计算机模型：回顾与机遇

【摘要】在现代科学中，计算机模型经常被用来理解复杂的现象，一个蓬勃发展的统计社区围绕着分析它们而发展起来。这篇综述旨在关注日益流行的随机计算机模型——为从业者提供统计方法目录，为统计学家（无论是否熟悉确定性计算机模型）提供介绍性观点，并强调与以下方面相关的开放性问题从业者和统计学家。高斯过程代理模型在本次调研中占据中心位置，并解释了这些模型以及随机设置所需的几个扩展。讨论中突出了设计随机计算机实验和校准随机计算机模型的基本问题。带有数据和代码的指导性示例用于描述各种方法的实现和结果。

【原文】 Baker, E. et al. (2022) ‘Analyzing Stochastic Computer Models: A Review with Opportunities’, Statistical Science, 37(1). Available at: https://doi.org/10.1214/21-STS822.

1 简介

计算机模型，也称为模拟器，到处都在使用。这些是描述和近似感兴趣过程的程序。该代码通常采用一组输入并产生一些输出。随机模拟器与确定性模拟器不同，由于存在随机元素，因此可以在相同输入的情况下产生不同的输出。这种计算机模型得到了广泛使用。例如，基于代理的模型 (ABM) 处理大量个体，其中在任何一个时间步长采取的特定行动具有复杂性和不确定性，不允许确定性建模。 ABM 很普遍（Johnson，2010；Johnson 和 Briggs，2011；Ramsey 和 Efford，2010；Smieszek 等，2011；Grimm 等，2006），用于探索社会学、交通、生态学、流行病学和其他现象

下面是一个随机模拟器实验的基本模型。如果代码在（向量）输入 $x$ 上运行，产生（标量）输出 $y(x)$ ，这可以表示为：

y(x) = M (x) + v, v\sim N (0, \sigma^2_v( x)), \tag{1.1}

其中 $M(x)$ 是输出的期望值 $E[y(x)]$ 。可变性 $v$ 说明了随机模拟器的随机性，最终由代码中的伪随机数生成引起。它的方差 $\sigma^2_v$ 通常取决于 $x$ ，常量方差作为特例包含在内。对于确定性模拟器， $\sigma^2_v = 0$ 。

随机模拟器中的随机性总是需要多次模拟，从而限制了可以有效处理的复杂性（包括输入维度的大小）。在随机模拟器中重复运行的前景引入了复制和探索之间的权衡，这是一个具有挑战性的设计问题。当噪声 $v$ 的方差不恒定时，噪声 $v$ 会对分析提出额外要求。本文研究了这些基本问题，确定了可访问且有效的方法，并指出了应该解决的未解决问题。

式（1.1） 通常用于模拟物理实验，其中观测值 $y(x)$ 为真值 $M(x)$ 加上测量误差（可能还有内在变异性），或者对于观测性研究，其中 $M(x)$ 适合具有残差的观测值。因为它们在结构上相同，物理实验可以用随机模拟器的方法进行分析（Gao 等, 1996）。然而，背景和目标往往不同，导致不同的问题表述和不同的结果解释。

方法的选择及其假设和局限性对于任何实验分析都是至关重要的。对软件简单性和可用性的倾向会鼓励对具有常数 $\sigma^2_v$ 的 $M$ 使用标准统计回归模型（例如，线性回归）。这种方法在某些情况下是有效的，尤其是当空间 $X$ 为可能的输入很小，引出了一个问题，即它作为一般处方的可靠性如何。由模拟器建模的复杂系统可能既不建议也不允许太多简化。本综述中描述的方法允许模拟数据在一般条件下指导方法的选择，几乎没有或没有简化。统计学（Sacks 等，1989；Kennedy 和 O’Hagan，2001）和应用数学（Sullivan，2015）在确定性计算机实验的设计和分析中发挥着重要作用。毫不奇怪，为确定性模拟器开发的一些方法进行了修改，可以在随机环境中使用。在许多情况下，由随机性驱动的备选方案是必要的。这些结构差异将在下面的叙述中注明。

1.1 目标

我们有三个主要目标；所有这些都与该主题的跨学科性质有关。

【目标一】：让学科科学家注意到有效的统计方法，并加深对使用中的随机模拟器的理解。 下面对使用（或引用）的统计工具的描述试图避免陷入数学复杂性的泥潭。包括个别方法的一些细节，以帮助理解这些方法的优点和缺点。在测试平台案例（ 第 2 节 ）中举例说明了多种方法的应用，并在可能的情况下确定了方法可用的软件。

【目标二】：是让统计学家熟悉一个对循证策略的形成至关重要的领域。 一般而言，基于代理的模型 (ABM) 和随机模拟器的研究和应用非常需要统计学家。熟悉确定性模拟技术的研究人员将看到直接的机会，但各种统计专业知识对于推进随机模拟器的研究至关重要。

【目标三】：推动统计学家、学科科学家和计算机科学家之间的密切合作。随机模拟器的分析是一个发展中的领域，有许多未解决的问题。挑战通常是由问题的规模和一系列问题驱动的，这些问题的解决需要统计学家、学科科学家和计算机科学家之间的密切合作。

本文结构如下：
第 3 节 描述了构成分析基础的模型；
第 4 节 专门讨论运行什么模拟器这一基本问题。
第 5 节 阐述了模拟实验的共同目标：校准。
第 6 节 讨论了其他重要的模型和目标，但更多的是在本次调研的 “边界” 上，因此不太详细。最后，
第 7 节 总结了结论并提出了未回答的问题。

本文的参考资料并未涵盖随机模拟器的全部工作，但连同此概述，应充分涵盖所讨论的问题。

2 示例的模拟器

本次调研将讨论三个随机模拟器以帮助理解。

其中两个被有意简化并用于展示所提出方法的关键特征。在某些情况下，更简单的策略可能同样有效，因为模型的复杂性已大大降低。由于使用的数据/生成机制可用，其他人可以比较不同的策略，但演示目的是与讨论和报告的计算相关的目的。

第三个是我们用来锚定和激励方法的模型。所讨论的具体模型是为应对 $2014$ 埃博拉疫情而开发的流行病学模型。对于埃博拉模型，合成的总体代表利比里亚的个体（总体约 450 万）及其活动时间表，引入时变接触网络个人和地点的概念，被开发出来（Mortveit 等，2015），并与基于代理的模型配对（Bisset 等，2009）。这个 ABM 模拟了在利比里亚从一个人传播到另一个人的传染病。由于传染性参数、可传播性仅在相互作用发生时控制感染概率，因此该模型（以及许多类似模型）具有随机性。该模型每天更新，疾病的进展由活动时间表、联系方式和其他流行病学特征决定。这个模型很复杂，具有高维输出、多个未知输入和非正态性。 Fadikar 等 (2018) 进行的分析使用了本文中讨论的想法，并解决了所有这些问题（参见 第 3.3.1 节、第 3.4 节 和第 5 节）。

2.1 鱼类捕获-再捕获

我们考虑的第一个简化的随机模拟器模拟鱼在标记-重捕获应用中的运动和鱼群行为。标记和重捕获涉及：捕获总体样本、标记样本、释放样本、再捕获样本（通过捕获另一批样本并计算有多少标记）。重捕获的数量被用于估计总体规模（Begon 等，1979）。该过程通过在具有边界条件的二维矩形湖泊中的随机位置初始化鱼群来建模。鱼开始根据简单的、基于代理的规则移动和学习。最初一段时间后， $100$ 条鱼在穿过湖中的 “网” 时被标记。第二阶段时间后，使用同一张网捕获 $100$ 条鱼，并记录其中 “重捕获” 的数量。

这种基于代理的模型是 NetLogo (Wilensky, 1999) 中开发的蜂群模型的修改版本。蜂群模型中出现的集体行为是为每个个体代理人提供同一组简单规则的结果 (Reynolds, 1987)。蜂群模型被修改为包括上述标记-重捕获动态。鉴于观测到的重捕获的鱼类数量，该模型可用于估计鱼类种群的总规模（见 第 5.3 节 ）。唯一考虑的输入是总种群中鱼的数量，输出是重捕获的鱼的数量。该模型的其他输入控制着鱼的个体运动规则；

为简单起见，这些在这里被忽略并设置为默认值。补充代码，以及编译好的Rmarkdown文档，对应我们对这个模拟器的分析，可以在 https://github.com/jhuang672/fish 找到。重新运行模拟器需要从 https://ccl.northwestern.edu/net\logo/ 安装 NetLogo。

2.2 海洋环流

第二个简化示例是一个随机模拟器，它模拟南大西洋薄水层（约 $2000$ 米深）中的氧气浓度（McKeague 等，2005；Herbei 和 Berliner，2014）。物理模型通过平流扩散方程（McKeague 等 (2005) 的方程 (4) ），即非线性偏微分方程 (PDE)，根据水流速度描述氧浓度的动态和扩散系数。

对于一组给定输入，平流扩散方程的解在封闭形式下不可用。不过，使用理论（Feynman，1948；Kac，1949）可以通过相关的随机过程（Herbei 和 Berliner，2014）得到近似解。对于域内的特定位置，会生成过程的随机路径，从而产生与该位置的 PDE 近似解的噪声结果。

通过使氧气浓度输出仅取决于四个输入来简化此示例：两个未知扩散常数（ $K_x$ 和 $K_y$ ）和两个位置变量（纬度和经度）。所有其他输入都固定在标称值。这种随机近似在物理科学中有很多，要么是由于计算限制，缺乏对底层系统的完全理解，要么是因为所研究的系统本身被认为是随机的。可以在 https://github.com/Demiperimetre/Ocean 找到与我们对该模拟器的分析相对应的补充代码和编译的 Rmarkdown 文档。

3 统计模型

运行模拟器并生成由式（1.1） 描述的数据输出的实验可以有多个目标。

一个主要目标，也是我们在这里关注的目标，是使用模拟数据来预测模拟器的值 $M (x) + v$ ，以及预测的不确定性，在未观测的 $x$ 处运行模拟器并不是免费的。当 $M$ 被认为是 “简单的” （例如， $x$ 坐标的多项式函数）时，可以使用许多标准的 “经典” 技术来近似 $M$ 。例如，线性回归模型和广义线性模型已被 Andrianakis 等 (2017) 使用和 Marrel 等 (2012)。诸如 第 2 节 中的那些复杂问题不太容易管理：为复数 $M$ 指定函数形式需要足够的先验知识或大量数据，而这两者通常都缺乏。本文的主要重点是为应对此类问题而开发的方法；对于更简单的情况，可以使用各种标准方法的充分参考。

在选择统计模型（这里称为代理模型，因为它充当计算机模型的代理）之前，需要考虑一系列因素。除了方法论假设之外，重要的是要考虑 “背景” ，即正在研究的特定问题的条件，从而得出 式（1.1） 及其扩展。一些重要的背景包括：

输入空间通常是一个超矩形：输入 $x$ 的每个坐标都受到上下限的约束。 第 2.2 节 通过采用矩形输入空间简化了问题，即使大西洋不是矩形。
式（1.1） 中的输出 $y$ 是标量，但多输出（例如时间序列）也很常见。
一些输入可能是分类的而不是数字的。
$v$ 的概率分布，变异性，通常被认为是正态的，但通常是无效的，就像埃博拉模型一样。

追溯到 Sacks 等 (1989); Currin 等 (1991) 的大量文献，主要是关于确定性模拟器的，发现高斯过程 (GP) 模型可以为 $M$ 生成一个灵活、有效的代理。

这种方法及其改进需要解决 式（1.1） 中依赖于输入的 $\sigma^2_v(x)$ ，可以有效地用于随机模拟，正如文献中记录的那样（例如 Kleijnen，2009，2017），并将在下面显而易见。在 O’Hagan (2006) 中可以找到全面的直观解释（对于确定性计算机模型）。从统计角度对高斯过程的更多技术描述可以在 Santner 等 (2018) 中找到和 Gramacy (2020)；有关机器学习的观点，请参阅 Rasmussen 和 Williams (2006)。

简而言之，GP 的使用允许计算机模型运行在选择替代品和评估其预测不确定性方面发挥关键作用。神经网络等深度学习方法和其他通用预测器也得到广泛使用。这些现代学习机器在产生不确定性和识别关键输入方面存在困难，但有针对此目的的积极研究（Neal，1996；Graves，2011；Welling 和 Teh，2011；Papamakarios 等，2019；Gal 和 Ghahramani，2016； Lakshminarayanan 等，2017）。

3.1 高斯过程代理

假设输入空间 $X$ 是一个d维的超矩形；输出 $y(x)$ 是单变量（标量）；并且这种变异性呈正态分布。此外，假设：

-A1：可变性 $v$ 具有常数方差 $\sigma^2_v$
-A2：平均值 $M (x) = \mu + Z(x)$

A3： $\mu$ 是常数
A4 $Z(·)$ 是 $X$ 上的高斯过程，均值为 $0$ ，协方差函数为 $K$ ，解构为方差 $\sigma^2_Z$ 和相关函数 $C$ 的乘积。

GP 的技术定义（ 假设 A4 ）是：对于任何有限 $N$ 和输入集合 $XN = (x_1,\ldots , x_N )$ , $ZN = (Z(x_1),\ldots , Z(x_N ))^{\top}$ 是均值为 $0$ 和 $N × N$ 协方差矩阵 $K_N$ 的多元正态随机变量，其条目为 $K(x_i, x_j)$ 。由此可见，模拟器输出 $Y_N = (y(x_1),\ldots , y(x_N ))^{\top}$ 也是多元正态分布，但具有均值 $\mu \mathbf{1}$ 和协方差矩阵 $K_N + \sigma^2_v\mathbf{I}_N$ ，其中 $\mathbf{I}_N$ 是恒等式 $N × N$ 矩阵， $\mathbf{1}$ 是元素值为 1 的 $N$ 维向量。

一种解释是，这些假设描述了均值 $M$ 的所有可能函数的先验分布。 GP 的不同选择允许不同类别的可能 $M$ ； GP 的强大之处在于这些类可以足够大以允许所有合理的可能性。在指定 $\mu$ 、 $K$ 和 $\sigma^2_v$ 之后，可以进行贝叶斯分析，从而得到所有函数的后验分布，这些函数在考虑到观测到的模拟器运行后仍然可以表示 $M$

$M$ 和 假设 A2 和 假设 A4 的另一种解释是将 $M$ 视为随机函数，其中 $\mu$ 是回归函数（如线性回归）， $Z$ 的 GP 对与 $\mu$ 的偏差进行建模。这两个公式具有相同的数学结构，但有不同的解释。

给定观测到的模拟器数据 $\{XN , Y_N \}$ 的任何新运行的预测分布 $y(x_{new})$ 也是正态的，并且具有已知的分析形式。预测的均值 $\mu_N (x_{new})$ 和方差 $\sigma^2_N (x_{new})$ 为：

\mu_N (x_{new}) = \mu + k_N(x_{new})^{\top}(K_N + \sigma^2_v \mathbf{I}_N)^{-1}(Y_N - \mu \mathbf{1}) \tag{3.1}

\sigma^2_N (x_{new}) = \sigma^2_v + \sigma^2_Z - (k_N (x_{new})^{\top}(K_N + \sigma^2_v \mathbf{I}_N )^{-1} k_N (x_{new}), \tag{3.2}

其中 $k_N (x_{new})$ 表示 $N$ -向量 $(K(x_{new}, x_1),\ldots , K(x_{new}, x_N ))^{\top}$ 所需预测和观测数据之间的协方差。一旦指定了相关函数 $C$ ，就可以从数据中估计参数（ $\mu$ 、 $\sigma^2_Z$ 和 $\sigma^2_v$ ）。为了指定 $C$ ，这里可以采用确定性模拟器采用的方法：指定参数化族 $C_{\theta}$ 并使用数据来估计 $\theta$ ，从而使 $C$ 适合观测。 $C_{\theta}$ 的一个例子是平方指数相关函数族（也称为高斯核）：

C_{\theta}(x, w) = \exp \big \{ - \sum\limits^{d}_{j=1} (x_j - w_j)^2 \theta_j \} \tag{3.3}

此相关函数适用于逼近维度 $d$ 上非常平滑、无限可微的函数。替代相关函数存在并被使用；一种常用的替代方法是 Matérn 5/2 相关函数（Stein，2012），它适用于逼近不太平滑的函数（只有 2 个导数）。

通过选择系列 $C_{\theta}$ 和 假设 A1-A4 ，观测到的输出的似然是可用的，并且可以计算最大似然估计 (MLE) $\hat{\mu}$ 、 $\hat{\sigma}^2_v$ 、 $\hat{\sigma}^2_Z$ 和 $\theta$ 。此后， $\mu_N (x_{new})$ 和 $\sigma^2_N (x_{new})$ 将用于表示预测分布的均值和方差，即使在估计方程 式（3.1） 和 式（3.2） 中的参数时也是如此。计算机模型输出 $y(x_{new})$ 的预测概率分布为：

y(x_{new}) \sim N (\mu_N (x_{new}), \sigma^2_N (x_{new})) \tag{3.4}

插入估计参数而不考虑其不确定性会丢失对不确定性的正确评估。因此，以这种方式获得的预测方差 $\sigma^2_N (x_{new})$ 称为插件（或标称）预测方差。在许多情况下，使用完整的贝叶斯分析来估计参数的替代方案在计算上可能不切实际，但并非不可能（中间方案和近似值已被证明是有用的，例如 Spiller 等（2014））。

对于 式（3.3） 中的相关函数，以及 Matérn 5/2 等其他函数， $Z(x)$ 和 $Z(w)$ 之间的相关性仅取决于 $x - w$ ，即两个输入向量之间的差异。也就是说， $Z$ 被假定为固定 GP（因此， $y$ 也是）。对于在输入空间的一个区域与另一部分表现出明显不同行为的函数，平稳性是有问题的。 Gramacy 和 Lee (2008), Ba 等 (2012)，Kersaudy 等 (2015), 和 Chen 等 (2016) 等解决和讨论了这个问题。 第 6.1 节 讨论了一种解决方案。

尽管上面的数学表达式相当复杂，但由于有许多可用的包（例如：R 中的 DiceKriging（Roustant 等，2018）、后面提到的 hetGP 包（Binois 和 Gramacy，2018），以及来自 Python scikit-learn 的 GaussianProcessRegressor 函数（Pedregosa 等，2011））。一般来说，GP 是一种灵活的方法，用于估计模拟器输出的平均 $M (x)$ ，尽管缺乏先验知识。这在 图 1 和 图 2 的顶部面板中进行了说明，但我们首先介绍了一个重要的建模转折，以应对随机计算机模拟的一个共同特征。

3.2 异方差高斯过程代理

常数方差 假设 (A1) 简化了统计模型的构建，因为只需要估计一个内在方差参数 $\sigma^2_v$ 。当 $\sigma^2_v(x)$ 被认为在输入空间内变化时，必须做更多的事情。 Boukouvalas 等 (2014a) 将 $\sigma^2_v(x)$ 建模为简单函数 $h$ （例如，多项式）的 $\exp(h(x))$ ，这是对仅假设一个方差参数 $\sigma^2_v$ 的简单扩展。（指数变换确保方差的正性。）与预测均值的类似方法（在 第 3.1 节 中简要讨论）一样，不清楚 $h$ 的用途，而且它的简单性可能无法满足许多应用中的复杂性。

几位作者将 GP 用于 $\sigma^2_v$ （Goldberg 等，1997；Kersting 等，2007；Boukouvalas 和 Cornford，2009；Ankenman 等，2010；Binois 等，2018a）。困难在于，这样做直接取决于观测输入 $X_n$ 处的 $\sigma^2_v(x)$ ，但这些值并未被观测到。如果在输入 $x_i$ 处有足够多的重复模拟运行 $r_i$ ，则 $x_i$ 处的样本方差（ $s^2(x_i) = \frac{1}{r_i - 1} \sum\limits^{r_i}_{j=1} (y(x_{ij})- \bar{y_i})^2$ ， $i = 1\ldots n)$ ）可用于估计输入 $X_n$ 处的 $\sigma^2_v(x)$ 。然后可以使用 式（3.1） 和式（3.2） 来预测 $\sigma^2_v(x_{new})$ 。（使用样本方差的对数，然后对结果求幂可以避免方差的负预测。）但是这种称为随机克里金法（SK，Ankenman 等，2010）的方法受限于需要足够数量的重复在每个输入以及分别处理方差和均值过程可能效率低下。

这些限制可以通过将输入的内在方差 $(\sigma^2_v(x_1),\ldots , \sigma^2_v(x_n))$ 作为未知参数（也称为潜在变量）以与所有参数相同的方式进行估计来消除其他未知参数。Goldberg 等 (1997) 以完全贝叶斯但计算量大的方式这样做。降低这些成本的努力构成了 Kersting 等 (2007) 以及 Boukouvalas 和 Cornford (2009) 方法的本质。Binois 等 (2018a) 提出的最新变体与可访问的软件 hetGP（Binois 和 Gramacy，2018）一起解决了计算风险，并且是本综述中描述和使用的方法。

解决异方差 GP (hetGP) 计算障碍的技术细节包含三个要素。第一，hetGP 将对数方差建模为 GP 在潜在（隐藏）变量上的平均输出。第二，使用 Woodbury 矩阵恒等式 (Harville, 1998) 将计算从处理所有 $N$ 个观测值减少到只涉及 $n$ 个唯一输入的计算，将计算复杂度从 $\mathcal{O}(N^3)$ 降低到 $\mathcal{O}(n^3)$ ，特别是当有许多重复。第三，元素使用 MLE 来设置所有参数。

虽然 Binois 等 (2018a) 提供了完整的细节。上面描述的第一个元素的一些细节值得注意。对于 $n$ 个不同的输入， $\lambda(x) = \sigma^2_v (x)/\sigma^2_Z$ 和 $\Lambda_n = (\lambda(x_1),\ldots , \lambda(x_n))，\log \Lambda_n$ 被视为 GP 的预测平均值对于潜在（隐藏）变量， $\Delta n = (\delta_1,\ldots , \delta_n)$ 。为便于说明，假设 GP 具有 $0$ 均值（常数均值实际上是 hetGP 中的默认设置）并将 $\Delta n$ 的协方差函数设为 $\sigma^2_g(C_g + gR^{-1})$ ，其中 $g > 0$ ， $R = \textrm{diag} (r_1,\ldots ,r_n)$ ， $C_g$ 是参数为 $\theta_g$ 的相关函数。然后 $\log \Lambda_n = C_g(C_g + gR^{-1})^{-1}\Delta n$ 。这种潜在的 $\Delta n$ 方法有助于平滑估计 $\Lambda_n$ 并为 $\lambda(x)$ 提供固定的函数形式，但不包含由于预测中固有 $\sigma^2_v(x)$ 的估计而产生的不确定性。给定 $\Lambda_n$ ，Woodbury 恒等式 (Harville, 1998) 降低了 $Y_N$ 的可能性，即所有输入（包括重复）的输出取决于大小为 $n$ 的数量。然后可以以 $\mathcal{O} (n^3)$ 的成本计算未知参数的最大似然估计。导数也可以以 $\mathcal{O}(n^3)$ 的成本计算，进一步促进了最大化可能性的优化。

附带说明一下，空间统计模型（通常与代理模型相关）中有时会出现异方差测量误差；然而，我们知道没有这样的模型允许以与 hetGP 相同的方式对内在方差过程进行完整建模和预测。例如，Nguyen 等 (2017) 的模型允许不同地点的非恒定测量误差，但它没有与其他模型参数一起估计这些测量误差，也不允许预测新的看不见的地点的测量误差。这主要是因为人们对预测空间统计中的测量误差过程兴趣不大（目标是 “真实” 的潜在信号），而对于随机模拟器，内在可变性可能具有直接的建模兴趣。

鱼的例子 。

我们将普通同方差 GP (homGP) 和 hetGP 代理应用于 第 2.1 节 中的鱼示例。模拟预算限制为 $400$ 次运行，重点关注种群中鱼类总数 $x$ 与第二轮捕获中重捕获的鱼类数量 $y(x)$ 之间的关系。总体是 $150$ 到 $4000$ 之间的整数。模拟器在 $[150,4000]$ 中的 $20$ 个唯一 $x$ 位置中的每个位置运行 $20$ 次，通过最大最小拉丁超立方体设计选择（参见 第 4 节 ）。计数的鱼数不能小于零，但正态性假设允许负鱼数，因此我们在执行分析之前对模拟输出进行平方根，然后对结果预测进行平方以返回到原始规模。此外，我们通过生成另一个数据集来估计 “真相” ；在相同的 $20$ 个站点中的每个站点复制 500 次。

应用具有平方指数相关函数的 homGP 代理产生 图 1 左上面板中的结果；右上面板显示了 hetGP 的结果。鱼模型的预测区间是在变换后的（平方根）空间中获得的，并平方以返回到原始空间。3 下图是具有 “真实” $2.5\%$ 、 $50\%$ 和 $97.5\%$ 分位数的图叠加。

关键结论是 homGP 和 hetGP 都捕获了非线性趋势（尽管在 $800$ 附近的区域有点偏离）。从真值图中可以清楚地看出非常数内在可变性的存在， $800$ 附近的区域显示出比其他地方更高的可变性。 hetGP 替代物并没有完全解决非恒定预测变异性，其中包括内在变异性和来自替代物的变异性，但确实改进了 homGP。全分辨率主要是模拟预算的问题，尽管替代设计可能会进一步改进 hetGP。我们的补充材料包括使用 第 4.3 节 的顺序设计方案改进的结果。外卖信息是 homGP 和 hetGP 都很容易处理这种趋势；异方差性鼓励使用 hetGP，可能需要添加模拟或改进设计。

海洋的例子 。

对于海洋模型（ 第 2.2 节 ），我们将每次模拟运行作为 $6$ 次模拟运行的平均值。已知真正的模拟器是非正常的；此调整使示例更符合高斯分布。现在，我们固定两个扩散系数， $K_x = 700$ 和 $K_y = 200$ ，将两个空间坐标作为唯一的变化输入。使用 $1000$ 次模拟（ $50$ 个站点，每个站点重复 $20$ 次），对于代理 homGP 和 hetGP，我们获得了预测平均表面和预测标准偏差表面（即，用于预测模拟器输出的标准偏差，同时考虑了不确定性围绕预测均值和内在方差估计 $\sigma^2_v$ ) 。这些表面绘制在 图 2 中，左列为 homGP，右列为 hetGP。

两个代理的平均曲面相似。 homGP 的预测标准差（左下）在整个输入区域相对恒定（明显受到内在方差 $\sigma^2_v$ 恒定的约束的影响）。 hetGP 的标准偏差表面明显不同，证明内在方差是非常数。 “真相” 是使用模拟器在 $500$ 个站点（通过 LHD 选择， 第 4 节 ）的重复运行（最多 $100,000$ 次）获得的，对每个站点的重复进行平均以获得真实平均值和与平均值的平方偏差以获得真正的方差）。这些绘制在附录（A 部分）和补充材料中；他们证实了非常数内在方差的存在。此外，hetGP 的标准差图显示出与真值图相似的结构，得出 hetGP 是该问题的更好替代品的结论。然而，这个结论有一个警告：由于设计和模拟的可变性（在 第 4.3 节 中进一步讨论），重复这个实验揭示了标准偏差图中的大量可变性。

总的来说，实现了对平均值的可靠预测，但不确定性不太确定。这类似于 Fish 示例，改善不确定性将需要更多模拟。这些结果表明 hetGP 优于 homGP。这通过 第 4.3 节 中的数值比较得到证实，其中还检查和比较了顺序设计。

3.3 非高斯可变性

在许多应用中，假设变异性 $v$ 服从正态分布是不合适的。例如，埃博拉病毒或鱼类模型中的计数数据是非正态的，并且不能小于 $0$ 。此外，给定输入 $x$ 处的 $v$ 分布可能不是单峰的：在埃博拉病毒示例中，即使输入 $x$ 固定的、重复的模拟可能导致两组不同的可能感染计数，这意味着双峰性。在某些模拟器中， $v$ 的分布可能更倾向于极值（肥尾）。有了这些可能性，正态性可能是一个需要谨慎使用的强有力的假设。

数据转换是一种历史悠久的方法，有时会在数据中产生 “足够” 的正态性，以允许使用基于高斯的方法（如 第 3.2 节 中的鱼类模型）。例如， Henderson 等 (2009) 使用 logit 变换 ( $\log y/(1 - y)$ ) 分析线粒体 DNA 中缺失的比例。 Plumlee 和 Tuo (2014) 通过关注输出分布的分位数采取了不同的路线——不需要正态性。这两种方法都具有导致对第 3.1 节 和 第 3.2 节 中的方法进行相对简单修改的吸引力。

还有更复杂的方法，通常缺乏同样的实施便利性。例如，Moutoussamy 等 (2015) 尝试对潜在的概率密度函数本身进行建模，而不是对输出 $y$ 进行建模。 Xie 和 Chen (2017) 设计了一个学生 $t$ 过程，它与 GP 过程没有太大区别，同时允许在数据分布中出现更重的尾部。

3.3.1 分位数克里金法

分位数克里金法 (QK) 是一种越来越流行的随机计算机模型仿真工具（Rannou 等，2002；Plumlee 和 Tuo，2014；Zhang 和 Xie，2017；Fadikar 等，2018）。这些方法是环境应用中使用的空间克里金公式（Zhang 等，2008；Zhou 等，2012；Opitz 等，2018）的自然延伸，通常会进一步定制模型以考虑罕见事件和极端分位数.

QK 方法直接对感兴趣的特定分位数建模，例如每个输入的中位数和下/上 $95\%$ 分位数。需要关于模拟器输出分布的最小假设。 $Q_q(x)$ ，输入 $x$ 处模拟器输出的第 $q$ 个分位数，用 GP 建模。给定输入 $x_1, \ldots, x_n$ 处的值 $Q_q(x_i)，\ldots$ ， $x_{new}$ 的分位数 $Q_q(x_{new})$ 可以使用 式（3.1） 和 式（3.2） 进行预测。该框架允许可变性 $v$ 的分布呈现几乎任何形状。虽然输出 $y$ 的真实生成过程丢失了，但我们可以描述它的分布。

要实施 QK，需要输入的目标分位数的值。正如在 第 3.2 节 中，可以使用输入的样本方差估计，这里可以使用样本分位数。假设每个 $x_i$ 处有足够的重复 $r_i$ ，则可以计算所述样本分位数。用于预测新分位数 $Q_q(x_{new})$ 的 GP 还应包括噪声项 $\sigma^_q$ ，以确认样本分位数是估计值。假设样本分位数的变异性呈正态分布也可能是无效的，但与感兴趣的数量 $y$ 相距更远，并且在实践中通常是可以接受的。

将分位数 $q$ 作为 GP 模型的附加输入包括在内可能是一个有用的修改。分位数 $Q_q(x)$ 可以重新表述为 $Q(x, q)$ ，将输入的维数从 $d$ 增加到 $d + 1$ 。该策略允许预测任何所需分位数 $q$ 的 $Q(x, q)$ ，而不是只是那些根据经验估计的，并被 Fadikar 等 (2018) 用于埃博拉模型。

基于 QK 的替代方法也在开发中。例如;一种称为非对称克里金法（AK、Zhang 和 Xie，2017）的有前途的 QK 变体通过利用分位数回归方法（Koenker 和 Bassett Jr，1978）不需要样本分位数。

鱼的例子 。

对于鱼模拟器，QK 是使用与以前相同的模拟数据集实现的，因为有许多重复可用。以20个总体规模的样本 $5\%$ 、 $27.5\%$ 、 $50\%$ 、 $72.5\%$ 和 $95\%$ 分位数组成观测数据，采用分位数q作为增加的输入维度进行修正。 图 1 显示了 $5$ 个不同分位数的预测 $Q(x, q)$ 均值以及数据（左图），并将 “真实” 值与 $5\%$ 、 $50\%$ 和 $95\%$ 分位数的预测值进行了比较（右图） ).

图 3 中中间的紫色曲线是预测的中值。外面的红线是预测的 $5\%$ 和 $95\%$ 分位数；内部蓝色曲线是预测的 $25\%$ 和 $75\%$ 分位数。 $5\%$ 和 $95\%$ 分位数的非单调 “波浪” 线反映了仅基于 $20$ 个观测值的极端分位数的自然变异性。如果没有大量的重复，准确捕获极端分位数是困难的，这是 QK 的一个缺点。呈现的其他分位数显示出更多的规律性。

QK 的结果与 图 1 中的结果差别不大，其中平方根变换就足够了。对于更复杂的问题，例如埃博拉模型，该方法适用，而其他方法可能不太适用。无论如何，如果有足够的数据来估计分位数，QK 可能是一个很好的稳健选择。

3.4 多输出

到目前为止的讨论都假定模拟器输出感兴趣的单个标量。对于多变量输出，最好使用更全面的模型。已经采用了构建多元 GP 的复杂方法（Conti 和 O’Hagan，2010；Fricker 等，2013；Paulo 等，2012）。还接受了针对时间序列和其他输出的定制的、特定于问题的公式（Farah 等，2014；Sun 等，2019）。如果输出数量较少，使用自己的代理模型独立处理每个输出通常就足够了。尽管忽略了不同输出之间的任何相关性并因此浪费了信息，但这种方法可能是有效的。例如，Spiller 等 (2014) 在一个地区的众多地点中的每一个地点部署独立的代理人，效果良好。

或者，通过将 $T$ 输出的索引 $t$ 视为附加输入维度（将输入空间的维度从 $d$ 更改为 $d + 1$ ），可以形成 $d + 1$ 维度上的 GP 代理（Bayarri 等， 2009).该方法允许对不同输出之间的相关结构进行建模。这类似于 QK 修改，其中分位数级别被视为添加的输入（ 第 3.3.1 节 ）。这种技术的一个缺点是，如果 T 非常大，则会出现计算问题，因为 GP 必须在 $N T$ 个数据点上进行训练，而不仅仅是 $N$ 。内在可变性阻止了 Bernardo 等 (1992) 使用的那种简化用于此设置中的确定性模拟器。

一种不同的方法通过使用基函数 $\psi(t)$ 来表示输出，从而将 $T$ 大小的影响降低到较小的 $K_0$ ：

y(x, t) = K_0 ∑ k=1 w_k(x)\psi_k(t ) + \delta(x, t)。 \tag{3.5}

系数 $w_k(x_i) k = 1,\ldots , K_0$ 由数据决定； $\delta(x, t)$ 是基函数表示和数据 $y$ 之间的残差。如果 $K_0 = T$ ，则 $\delta = 0$ 。通常， $K_0$ 被认为远小于 $T$ 但足够大，因此误差 $\delta$ 足够小。每个 $w_k(x)$ 都可以用一个替代项独立建模， $y(x, t)$ 的预测从 式（3.5） 中获得，忽略 $\delta$ 。

不同的基选择可能适用于不同的设置。例如，Bayarri 等 (2007a) 在 $t$ 是时间的确定性设置中对 $\psi_k$ 使用小波。基函数的一个常见选择是主成分： $\psi$ 是矩阵 $Y ^{\top}_N Y_N$ 的特征向量，其中前 $K_0$ 个特征值与前 $K_0$ 个特征值按降序排列。通常情况下，前几个（五个或更少）主成分足以捕获有关完整 ( $T$ ) 数据集的足够信息。系数 $w_k(x_i)$ 则等于 $\sum^{T}_{t=1} y(x_i, t) \psi_k(t)$ 。有关主成分的更多信息可以在 Jolliffe (2011) 中找到，用于获取 $\psi$ 和 $w_k$ 的软件很普遍。

关于使用主成分对高维模拟器输出建模的进一步讨论可以在 Higdon 等中找到。 (2008)。 Fadikar 等也使用了主成分。 (2018) 对随机埃博拉模拟器的时间序列输出进行建模。虽然主成分是常见的默认值，但人们担心数据集的关键特征可能会留在丢弃的 $\delta(x, t)$ 中，从而妨碍可靠的预测。Salter 等 (2019) 记录了这些关于校准的担忧并提出了替代方案。

对于函数输出的问题，可能丢失数据和/或不规则间隔的数据（例如不规则间隔的时间步长或空间位置），函数分解也很有用。例如，Ma 等 (2019) 使用功能主成分分析来模拟卫星观测模拟。

4 实验设计

对于实验，设计（x 值的选择）和分析（输出 y(x) 的评估）原则上是密切相关的。其他的考虑也可以进入。对于物理实验，通过分块和随机化来控制外部影响或干扰因素通常是设计的重要组成部分。计算机实验中不存在外部影响，因此控制有害因素通常是无关紧要的。然而，许多次要参数通常是固定的，可以改为随机化，从而增加了固有误差。

具有特定目标（例如，预测模拟器输出）和准确性标准（例如，平均预测不确定性：集成均方预测误差，IMSPE5），优化标准的设计是首选。由于标准通常取决于代理，而代理又取决于未知参数，因此在收集任何数据之前使用什么作为参数的替代是一个问题。对单阶段确定性计算机实验的广泛研究通过淡化最优性并推荐易于计算的 “空间填充” 设计解决了这一难题，这些设计不会遗漏大面积的输入空间。空间填充设计很容易计算，而优化 IMSPE 很复杂且没有实质性优势。为了实际采用，设计必须易于生产且有效。

存在多种获得空间填充设计的方法，最流行的是拉丁超立方体设计（LHDs；McKay 等，1979）。6 LHDs 已被证明是足够的，特别是当与附加标准结合使用时，例如最大最小标准，其中一个还最大化了设计中点之间的最小距离。7 即使是随机的 LHD 通常也足够了。 Sobol 序列设计 (Sobol, 1967) 对于预测确定性模拟器的输出同样有效。 Sobol 设计的 xs 是按顺序生成的，这使得在使用多阶段或顺序设计策略时很容易保留空间填充特性。8 Pronzato 和 M uller (2011) 对这些和其他空间填充进行了冗长的讨论方法，一些与非矩形几何有关。

对于随机模拟器，情况远没有那么清晰。内在可变性的存在增加了复制的复杂性，这在确定性实验中是不存在的。使用相同的输入，随机模拟器可以多次运行（复制），由于固有的随机性，每次都提供不同的输出值。重复显然对内在方差 \sigma^2_v 的估计有影响，因此对预测有影响（见 第 3.2 节 ），因此重复的数量和位置很重要。单阶段实验的一种简单方法是使用空间填充设计来建立实验的站点 X_n = (x_1,\ldots , x_n)，然后在每个站点添加重复。确定数量，ri，每个站点 x_i 的复制以及如何在复制和站点之间分配，即如何选择唯一站点的数量，n，给定总模拟预算 N，不是很好理解。事实上，尽管有数字证据并且人们普遍认为复制可能是有利的，至少在适当的情况下，但完全需要复制的理论证据有限。例如，Wang 和 Haaland (2019) 通过最小化 IMSPE 的边界来生成设计。他们的数值结果表明不需要重复，除非 \sigma^2_v(x) 与 \sigma^2_Z（估计均值 M 时测量不确定性的一个因素）相比很大。

内在可变性的存在表明在多阶段设计中存在价值，其中第 1 阶段用于获取有关 \sigma^2_v(x) 的信息，后面的阶段利用此信息来分配重复和选择新输入。问题是如何为第 1 阶段选择输入，以及如何利用第 1 阶段的结果来选择新的输入并在以后的阶段进行复制。复制和多阶段（包括完全连续的）这两个因素是制定适当设计策略的核心。在下面的讨论中注意这两个因素。

4.1 单级设计

单阶段研究中的一种常见方法是对输入使用空间填充设计，比如数量为 n，并且在每个输入处重复 r，有时没有重复，即 r = 1。根据 第 3 节 中的描述进行预测选择的特定预测模型。必须对运行总数和每个输入站点的重复次数进行选择 (N = nr)。通常，N 是一个预算问题，但很少有人深入了解应该如何选择 r，除非满足特定的代理模型要求，如 SK（ 第 3.2 节 ）。

对于他们的单阶段研究，Marrel 等 (2012) 使用没有重复的标准 LHD 来比较不同统计模型的性能。另一方面，Plumlee 和 Tuo (2014) 使用的 LHD 在每个 x_i 处具有不同数量的重复 ri。在他们的情况下，重复的数量必须很大，因为 QK 方法（ 第 3.3.1 节 ）取决于设计的每个输入站点的输出 y(x_i) 的计算分位数。

4.2 两阶段设计

两阶段设计的情况主要是为了能够在阶段 1 估计 \sigma^2_v 并将其用于第二阶段。安肯曼等 (2010) 在 SK 的背景下提供一种解决方案。第一阶段设计通过大小为 n1 的 LHD 选择 x_is，每个输入的重复次数为 r，从而导致第 1 阶段的运行总数为 N1 = n1r。第一阶段分析使用r 在每个输入处重复使用样本方差来估计 \sigma^2_v (x_i)。如 第 3.2 节 所述，然后使用 GP（使用 \log s2(x_i)）为所有 x 生成 \sigma^2_v(x) 的 “插件” 估计。不同的 GP 使用该方差估计来构建平均输出 M 的预测器。

对于阶段 2，选择了 n2 个额外的唯一输入位置，以便组合的设计位置集 X_n = (x_1,\ldots , x_n) 保持空间填充。然后使用阶段 1 中构建的 GP 模型，通过整合 MSPE 所有可能的输入 X 来计算 IMSPE。根据重复次数 Rn = (r1,\ldots , rn) 最小化 IMSPE 可提供最佳重复次数对于所选的 X_n。详情见 Ankenman 等（2010）。一个困难是，最佳 Rn 可能会为第一阶段站点生成一个 r_i，该 r_i 小于第 1 阶段已经使用的 r。然后需要对该方法进行一些修复。

在此设置中，Sobol 序列可用于获得在第 1 阶段和第 2 阶段都进行空间填充的设计。这不是 Ankenman 等所做的。 (2010)，但是 Sobol 序列更容易实现并且可能产生类似的结果。也可以通过优化 IMSPE 为阶段 2 选择唯一输入 X_n，但会增加计算负担。缺乏对 n1、n2、N 的值和每个不同输入的重复的适当建议（在 Ankenman 等 (2010) 中，建议是临时的），并且对于单阶段实验，开放研究。可以通过重复上述过程中的第 2 阶段来构造第三阶段设计（或者实际上，任何多阶段设计）。

4.3 时序设计

当统计设计和结果分析密切协调时，执行顺序过程可能是可行的，在第一阶段之后，一次选择一个运行。每次运行后，可以更新所有相关数量以确定下一次运行。这解决了在不预先指定其分配的情况下了解 \sigma^2_v 和获取新运行的问题。顺序设计的一个优点是可以在达到预算约束之前满足标准时停止。另一个优点是增加了对模拟器、复制或其他方式进行有用运行的可能性。对于某些目标，例如优化（ 第 6.3 节 ），顺序设计通常是必不可少的。对于全球预测，Binois 等 (2018b) 提出了一种顺序设计方法，在前面提到的 hetGP 包中实现。

Binois 等的策略。 (2018b) 从第 1 阶段开始，采用 n1 输入的空间填充设计 D1 和运行分配 (r(x_1),\ldots , r(x_n1))。如 第 3.2 节 所述，使用 M 的 GP 和先验 \sigma^2_v 的潜在 GP，MLE 计算处理所有参数，导致预测变量和 IMSPE(D1) 的可计算估计。考虑新点 z，作为新的唯一输入 x_n1+1 或作为 D1 中现有输入的复制。如果 z 最小化 IMSPE(D1 + z)，则将选择 z 添加到设计 D1，从而产生新的设计 D2。可以迭代此近视规则，每次添加新点时，代理项（包括其参数的 MLE）都会更新。当满足标准或计算预算耗尽时，该过程停止。

计算可行性因每次运行后所需的更新而紧张。另一方面，它的 “贪心” 性质减轻了计算负担：它只为下一次模拟器运行寻找最佳数据点，而忽略从长远来看可能更好的运行。

这不是唯一可用于全局预测问题的顺序设计方案。例如，TGP 和 BART（参见 第 6.1 节 ）中使用的树生成过程提供专门的顺序设计策略。 Gramacy 和 Lee (2009) 以及 Chipman 等提供了详细信息。（2010）。

模糊了多阶段设计和顺序设计之间的界限，有时分批运行额外的模拟是可行的（例如，在高效使用多核超级计算机时）。在这种情况下， “批量设计” 是可取的。这些是为确定性实验开发的（Loeppky 等，2009a；Duan 等，2017；Erickson 等，2018），但尚未明确扩展到随机情况。

当完全顺序的方法可行时，上面概述的 seqhetGP 策略是有价值的。有几个方面值得考察：

• 在设计的构建中广泛使用替代品需要通过评估替代品质量的诊断进行调研。 • 序贯策略的第一阶段必须避免糟糕的（例如，太小的）初始设计，以免糟糕的起始代理导致之后的糟糕选择。 • 顺序设计的效用取决于与仿真器运行相比实施的相对成本。对于具有挑战性的问题，模拟器运行的成本可能足以使顺序设计具有吸引力。 • 可以对顺序设计进行修改，以减少计算量，而无需付出大量准确性成本。例如，定期而不是在每一步之后重新估计参数。

海洋的例子。对于海洋模型，我们使用 50 个站点的初始设计，由 2-d 中的最大最小 LHD 选择，每个站点有 5 个重复。然后通过顺序方案分配剩余的 750 个数据点。生成的均值和标准偏差曲面如 图 4 所示。对于标准偏差曲面，设计站点与在站点进行的复制次数一起叠加。左侧面板中的平均表面与非顺序分析略有不同（ 图 2 ，顶行）。标准差表面看起来非常不同。对于设计本身，新输入在标准偏差较大的区域被大量复制，而在标准偏差较小的区域则较少。此外，顺序设计包括比固定设计更多的独特站点，以及输入空间边界上的更多点。

使用 第 3.2 节 中建立的 “真相” ，我们可以比较三种方法的性能。如前所述，异方差性的视觉存在阻碍了使用 homGP。在视觉上区分 hetGP 和 seqhetGP 代理的性能更加困难：均值看起来相似，虽然真实标准偏差中的某些模式似乎被 hetGP 捕获，但缺陷是可见的，并且幅度并不总是正确的。使用 seqhetGP 标准偏差，细微差别似乎丢失了。为了正确比较不同的方法，数值比较可能更有价值。

两个有用的数值度量是均方根误差、RMSE（代理人预测的均值与 “真实” 均值之间的平均平方差的平方根）和分数（Gneiting 和 Raftery (2007) 中方程 27 的正确评分规则） ). RMSE 衡量均值预测的准确性，而 Score 是测试组合均值和方差预测准确性的总体度量。使用一组测试输入 x_1,\ldots , xp 和模拟器输出 y1,\ldots , yp, 替代预测意味着 \mu_N (x_1),\ldots , \mu_N (xp) 和方差 \sigma^2_N (x_1),\ldots , \sigma^2_N (xp), 得分为 1 p p ∑ i=1  - ( y_i - \mu_N (x_i) √\sigma^2_N (x_i) )^2 - \log(\sigma^2_N (x_i))   。 \tag{4.1}

RMSE 越小越好，而 Score 越大越好。

对于这三种方法，homGP、hetGP 和 seqhetGP 的 RMSE 分别为 2.056、1.985 和 1.567；分数分别为-3.999、-3.880 和-3.834。 RMSE 结果表明 seqhetGP 最擅长预测均值，这在图中并不明显。 hetGP 和 seqhetGP 的分数接近但明显好于 homGP，证实了异方差性的存在。

随机模拟器中的随机性以及设计中的可变性（有许多可能的最大最小 LHD）会导致特定结果的很大程度的可变性，例如刚刚引用的结果。因此，很难依靠单一结果进行比较。因此，上述实验重复 100 次，结果 100 RMSE 和分数总结在 图 5 的箱线图中

箱线图证实了单一数据集的发现：存在异方差，seqhetGP 是首选。对重复实验的许多标准偏差图进行目视检查（如 第 3.2 节 中所讨论和在补充中发现的那样）揭示了相当大的变化和与真实标准偏差的偏离。如果没有丰富的数据，方差（内在方差和均值的 GP 不确定性）可能很难得到正确的结果，并且使用未考虑不确定性的插件估计会加剧困难。

4.4 统计模型参数估计设计

第 4.2 节 和 第 4.3 节 构建依赖于基于第 1 阶段数据的替代模型的设计，以便选择后续数据点。设计的质量取决于代理的准确性，而代理的准确性又取决于其参数的准确性。 4.2 和 第 4.3 节 中使用的方法的另一种方法是构建一个初始设计，其明确目的是更好地估计这些参数。

Boukouvalas 等 (2014a) 解决这个问题并专注于 hetGP 模型，使用一个简单的方差参数函数 (\sigma^2_v(x) = \exp(h(x))，其中 h 是一个简单的函数（例如，多项式）。他们提出最大化先前由 Abt 和 Welch (1998) 用于确定性模拟器的标准的设计：Fisher 信息矩阵行列式的对数，\log |I|。数值结果表明该方法在估计参数方面有所改进，但总体全局预测并不比使用空间填充设计更好，有时甚至更糟。当预测至关重要时，问题就出现了，如何在多阶段或顺序设置中将此类设计用于第一阶段，它对可以感觉到获得更好的初始代理模型。例如，参见 Zhang 等 (2020)。

5 校准

当模拟器的输入既未知又不可测量时，需要进行校准，这是实践中的常见情况。埃博拉模拟器中的传播率和海洋模型中的扩散系数就是此类输入的示例。为了（间接）推断这些输入的值并产生预测，需要以现场数据（实验或其他方式）形式添加信息。包含现场数据和校准参数，标记为 uC，导致观测模型：

yF (x) = yS(x, uC) + \deltaMD(x) + , \tag{5.1}

其中 yF (x) 是在可控（或可测量）输入 x 处的真实世界现场观测，yS 是具有额外未知、不可测量输入的模拟器 uC 是观测值 yF (x) 的测量误差（方差为 σ2） , 而 \deltaMD(x) 是一个重要的术语，它解释了模拟器不是现实的完美代表。 yF 错误地 “观测” 现实；现实 = yS + \deltaMD。

校准存在多种相互竞争的方法甚至哲学。下面概述了校准问题的几种解决方案。尽管校准在计算机实验中处于中心地位，但缺乏全面的比较。

5.1 Kennedy-O’Hagan 校准 (KOH)

方程式 5.1 中的公式是由 Kennedy 和 O’Hagan (2001) 为确定性模拟器制定的，是此后大部分校准和相关预测工作的基础。 Kennedy 和 O’Hagan (2001) 所追求的策略，在 Bayarri 等中实施。 (2007b)，获得了 yS 的替代项，并使用 GP 对 \deltaMD(x) 建模（尽管其他选择也是可能的）。用替代项替换 yS 后，可以通过贝叶斯分析获得所有未知数的后验分布。在实践中，替代模型仅使用模拟器数据进行拟合，忽略了现场数据可能产生的影响。这种模块化方法的详细信息和讨论可以在 Bayarri 等中找到。 (2007b) 和 Liu 等 (2009)

KOH 方法强调同时解决校准和模型差异的必要性。 uC 和 \deltaMD(x) 之间不可避免地会发生混淆，因为 uC 和 \deltaMD(x) 的多种组合会导致相同的实地观测数据。因此，uC 是不可识别的，其估计会受到影响，差异也是如此。尽管如此，对 y 和 E(y) 的最终预测是合理的，即使对 uC 和 \deltaMD(x) 的单独估计不是。有关详细信息和进一步讨论，请参阅 Higdon 等 (2004), Bayarri 等 (2007b)、Brynjarsd ́ ottir 和 O’Hagan (2014)，以及 Tuo 和 Wu (2016)。

已经出现了多种规避混杂的尝试。拓等 (2015) 通过正式将其定义为最小二乘量来减轻 uC 中的歧义； Gu 和 Wang (2018) 针对 Tuo 等之间妥协的差异提出了新颖的先验。 (2015) 战略与 KOH；和 Plumlee (2017) 介绍了与先验均值正交的差异先验。在随机模拟器文献中，Oakley 和 Youngman (2017) 移除了 \deltaMD，但通过夸大 uC 的先验分布中的可变性来进行补偿。完全忽略 \deltaMD 可以通过模拟器准确的有力证据来证明，但这种证据很少见。

对于随机问题，现实是随机的，差异项 \deltaMD(x) 不能假定为确定性的。对差异建模可能会受到模拟器模型的影响，同时认识到差异通常更平滑。例如，如果对 yS 建模需要具有 Matern 5/2 相关函数的 hetGP，则可能需要 hetGP 来解决差异，也许具有更平滑的平方指数相关性。在这种情况下进行完整的贝叶斯分析可能会非常昂贵，并且必须修改上述程序。宋等。 (2019) 使用 hetGP 来解决差异（但对于确定性模拟器），通过最大似然估计参数并遵循 Tuo 等 (2015) 以避免混淆。

重温埃博拉。埃博拉研究（Fadikar 等，2018）使用 KOH 框架校准了 ABM。模拟器 yS 有 5 个未知的、未测量的输入 uC，输出是截至第 1 周以及此后每周最多 57 周的受感染个体累计数量的日志。现场数据 yF 是一组报告的累积计数。对于统计模型，QK 策略（ 第 3.3.1 节 ）之后是将每个不同的模拟复制 100 次，然后在每个时间点（具体来说， $5\%$ 、 $27.5\%$ 、 $50\%$ 、 $72.5\%$ 和 $95\%$ ）压缩成均匀分布的分位数分位数）。然后使用 第 3.4 节 中概述的主成分分解将这些分位数输出轨迹减少到更易于管理的 5 维。

该方法的基础是假设流行病轨迹（实际和模拟）可以通过分位数轨迹来近似（即，在时间 1 类似于第 qth 分位数的已实现流行病在以后也类似于第 qth 分位数）。因此，分位数 q 作为输入参数（参见 第 3.3.1 节 ）包含在内，以允许 KOH 校准了解 5 个校准参数以及观测到的流行病的（未知）q 值。因为模拟器和现实分位数轨迹之间的差异不可能是嘈杂的，所以差异被视为确定性的（使用平滑样条而不是 GP 来处理差异）。获得未知 uC、\deltaMD 和 q 的后验分布，并用于预测累积计数和其他量。

在将现场数据限制为前 20 周的主要分析中，模型差异的估计几乎为零。使用第 42 周之前的现场数据进行的后续分析暴露了模拟器的一些不准确性（非零 \deltaMD(x)）——模拟器继续预测感染，即使在现实中流行病已经死亡之后。

海洋的例子。之前的海洋分析固定了两个扩散系数。实际上，它们是未知的，需要校准。 “现场” 数据是通过使用先前固定的扩散系数值（K_x = 700 和 K_y = 200）在 150 个不同的经纬度坐标下对 200 多次模拟进行平均而人工创建的。 “真” 值是通过添加假差异获得的，作为来自具有平方指数相关函数的 GP 的单一实现，方差为 1.64，\theta 值为 (1, 2)（ 式（3.3） ）。在这些基础上，添加方差为 4 的正态分布假设 “观测误差” ，每个站点有两个这样的观测结果。在实际问题中，现场数据将被观测而不是像这样生成。请注意，此问题的现场数据对应于模拟器的平均值，而不是来自模拟器的个别抽取；模拟器是随机近似的结果。

由于扩散系数现在不确定，模拟器有四个输入。计算机实验设计为在用于现场数据的 150 个站点和校准参数 K_x 和 K_y 的 500 个独特选择上运行。这是通过将 150 个经度和纬度站点的副本与大小为 500 的 max_imin LHD 相结合来完成的（ K_x，K_y），然后通过最大化4维空间中设计点之间的最小距离来改进组合设计。将这组点称为 Doc。通过在 Doc 中的每个点重复 10 次来执行模拟器实验。两个不同的代理（一个 homGP 和一个 hetGP）适合这个固定设计。此外，构建了一个 seqhetGP 代理，初始设计只有 4 个 Doc 副本，其余预算按照 Binois 等的策略分配。 (2018b)，如 第 4.3 节 所述。

对于以模块化方式完成的 KOH 分析，替代项仅适用于使用模拟数据。因为这里的现实是由模拟器的期望（而不是模拟器输出本身）表示的，所以 式（5.1） 中的 yS 被替换为 E(yS)。类似地，因为现实是确定性的，\deltaMD 被建模为标准 GP。当然，模拟数据是 yS 的输出，而不是 E(yS) 的输出——代理项用于近似确定性 E(yS)。然后使用 MCMC 获得剩余未知数的后验分布：扩散系数 K_x 和 K_y；模型差异GP、σ2 MD和\thetaMD的方差和相关参数；和观测误差 σ2)。关键参数的后验分布如 图 6 所示；它们的真实值为 K_x = 700、K_y = 200、σ2 MD = 1.64 和 σ2 = 4。

对于所有三个替代模型，K_x 的后验分布相当分散。 K_y 后验值高度集中，但不完全围绕真实值。观测误差的三个后验概率非常相似，但都指向更接近 5 而不是真实的 4 的估计值。差异方差的后验概率是弥散的。这些图强调了校准的困境：在存在模型差异的情况下获得准确的校准（和其他）参数值是有问题的。此外，对于嘈杂的数据，很难获得精确的估计。然而，KOH 确实产生了有用的后验预测分布。

表 1 将 KOH 校准的预测与其他 3 种校准方法进行了比较。第一个通过普通最小二乘法 (OLS) 估计 K_x 和 K_y：选择 (K_x, K_y) 使得平均代理预测和观测数据之间的残差平方和最小化。然后通过用参数 (K_x, K_y) 替换为 OLS 估计值 ( ^ K_x, ^ K_y) 运行代理来预测新的观测结果。第二种方法遵循一种经常采用的做法，即猜测或 “明智地选择” K_x 和 K_y 的特定值。这里，选择 K_x = 600 和 K_y = 400，然后使用代理项进行预测。将此方法称为 S\mathbf{I}_NGLE。第三种方法 NOCAL 生成预测，就好像没有现场数据一样，并且 (K_x, K_y) 的分布被视为它们的先验分布，[100, 1000] 上的独立均匀先验。在这些替代方法中，观测误差方差固定为真实值，并假设差异为 0（前者过于慷慨，后者在实践中太普遍了）。对于 NOCAL，氧气浓度的分布是通过从它们的先验分布中采样 K_x 和 K_y 的值并将它们插入替代项中获得的，而对于 KOH，通过从所有未知数的后验分布中采样来获得。

尽管 RMSE 的差异可以忽略不计，但分数表明 KOH 表现最好。 OLS、S\mathbf{I}_NGLE 和 NOCAL 的准确性也可能被夸大了，因为观测误差方差被认为是已知的，而在 KOH 中它是估计的。最小二乘估计值 ( ^ K_x, ^ K_y) 并不总是接近真实值，这并不奇怪，因为存在差异，以及替代项中可能存在的缺陷和高可变性。出于类似的原因，三个代理人之间几乎没有差异。

RMSE 的相似性是由于代理变量的巨大可变性、差异的存在、经度和纬度输入的主导以及校准输入的微弱影响的结果。第一个解释了 RMSE 的大小，最后一个解释了为什么相当不准确的校准输入（在 OLS 和 S\mathbf{I}_NGLE 中）无关紧要。因为 KOH 解决了差异，它的分数超过了其他的，表明对差异的解释是必要的，不能随心所欲

估计校准参数这里的重点与 第 4 节 一致，一直是改进全局预测。如果问题是在没有模型偏差的情况下为校准参数提供良好的估计，那么不同的设计可能更适合。丹布林等。 (2018) 在确定性模拟器的背景下解决了这个问题，但尚不清楚这些方法如何扩展到随机模拟器。此外，虽然 KOH 有助于进行有能力的预测，但 KOH 的复杂性和记录的缺陷导致了竞争性校准技术，这些技术也很常用。

5.2 历史匹配（HM）

历史匹配 (HM) 是 KOH 校准的常见替代方法（Craig 等，1997；Vernon 等，2010；Boukouvalas 等，2014b；Andrianakis 等，2017）。 HM 搜索模拟器输出与观测数据密切匹配的输入，同时识别各种不确定性的存在，包括模型差异。 HM 方法以直接的方式排除 “不可信” 的输入，而不是试图找到可能的输入。通过观测 yF ，并最初假设 uC 构成所有未指定的模拟器输入，如果满足以下条件，则 uC 被认为是不可信的： |yF - \mu_N (uC)| √\sigma^2_N (uC ) + σ2 MD + σ2 ) ≥ 3, \tag{5.2}

其中 \sigma^2_N、σ2 MD 和 σ2 分别是代理方差、模型差异和观测误差。换句话说，如果相对于这些不确定性，观测结果与使用该输入的模拟器输出之间的差异足够大，则该输入是不可信的。数字 3 来自 Pukelsheim (1994)，他表明至少 $95\%$ 的单峰分布包含在三个标准差内。当有多个输出或额外的可控输入时，方程 5.2 会有所修改（Vernon 等，2010）。

这个过程可以在所谓的 “波浪” 中重复，使用在一个波浪中发现的非难以置信的 uC 来为下一波浪生成模拟运行，依次减少 uC 可能存在的空间。通过这些波，HM 旨在避免 uC 不太可能出现的输入区域，在这方面，HM 是一种校准设计方案。在任何给定的波浪中，uC 的所有值都可能被认为不可信——所谓的终端情况（Salter 等，2019）——通常意味着 σ2 MD 设置得太低或模拟器不适合目的。安德里亚纳基斯等 (2015) 包含对 HM 的详尽描述，同时将其应用于复杂的 HIV 流行病学模型。

HM 和 KOH。对于 KOH，uC 的估计因差异而混淆，但预测及其不确定性是可用的。然而，在复杂问题中实施 KOH 即使不是棘手的，也可能是繁重的。推测上，一种混合策略可能是使用 HM 来减少输入空间，确认不存在终端案例，然后在缩小的空间中应用 KOH 来获得预测和不确定性。与本综述中的鱼类和海洋示例不同，复杂模型将是这种方法最有吸引力的模型。这种混合策略是进一步探索的主题。

5.3 近似贝叶斯计算（ABC）

获得校准参数 uC 和预测的后验分布（例如， 第 5.1 节 中的 KOH 方法）在计算上可能具有挑战性。 ABC 方法提供了一种替代方法，已被发现在中等复杂的环境中很有用（Rutter 等，2019）；但在更雄心勃勃的环境中则不然（McKinley 等，2018）。

ABC 的目标是从 π(\theta|YF ) 中生成样本，即未知数 \theta 的后验分布，给定现场数据 YF 。对于校准，将 \theta 视为 uC。 ABC 通过为未知数 \theta(s) 和来自 π(YF |\theta)π(\theta) 的输出 z(s) 生成样本来做到这一点，也就是说，从给定未知数的数据的可能性乘以先验概率的未知数。对于计算机模型，从似然生成样本相当于运行模拟器。此类样本仅在 z(s) = YF 时才被接受。对于无法出现完全相等的连续设置，如果 B(z(s), YF ) < τ 则改为接受，其中 B 是距离的度量，而 τ 是公差水平。然后通过接受的 \theta(s)s 的集合给出近似后验分布。当有多个输出时（或有其他可控输入 x，因此对于任何给定的 \theta(s)，实际上有多个输出），YF 和 z(s) 可以用信息汇总统计代替。找到足以满足所有输出的单一统计数据具有挑战性，而且选择不当的统计数据可能会使结果无效。

公差 τ 的选择很重要。如果 τ 很小，那么可能需要很长时间才能生成满足不等式的单个样本。如果 τ 不小，则对后验的近似不太可靠。对于校准，τ 可以解释为观测误差和模型差异的界限，从而导致 “正确” 的后验而不是近似（Wilkinson，2013）。这类似于 HM 与边界的主观选择。

ABC 可以在不使用代理的情况下完成，但可能需要多次运行模拟器本身来生成许多 \theta(s)。否则，可接受的 \theta 将保持太少，或者将需要 τ 的值过高。在任何一种情况下，准确性都会受到影响。这种计算障碍可以通过使用代理来缓解。

鱼的例子。在这里，我们将 ABC 应用于鱼类模拟器，以估计种群中有多少条鱼。假设在第二轮中捕获了 25 条鱼。确定总种群大小的一种直接方法是从 NetLogo 鱼类模型中多次模拟总鱼类种群的许多不同值，并在每次导致 25 条鱼被重捕获时 “接受” 模拟。这正是 ABC，并且是 ABM 的一种相当普遍的做法。这样做 10,000 次，对 200 到 4000 之间的整数使用统一先验，因此每个这样的总体规模都有 1/3801 的先验概率，产生 图 7 左侧面板中的结果。这种直接使用模拟器只产生 52 个接受样本，这是一个非常小的数字，这来自 10,000 个模拟器运行。相比之下，仅运行 400 次的 hetGP 代理拟合，可以从中快速抽取 1,000,000 个样本，产生 3811 个可接受的样本。这个结果如 图 7 的右面板所示，给出了具有相同整体形状的噪声较小的直方图。如果基于主体的模型成本甚微，那么代理对于 ABC 计算无疑是有价值的。

5.4 相关校准技术

Bound-to-Bound (Frenklach 等, 2016) 类似于 HM，其中扫除所有不确定性的误差边界被类似地定义，然后使用二次规划来寻找 uC 的可行边界。贝叶斯融合 (Poole and Raftery, 2000; Raftery 等, 1995) 是一种与贝叶斯校准相关的技术，用于调和模拟器输入和输出的引出先验分布之间的差异。它已应用于生态学、流行病学、城市建模和污染监测（ˇ Sevc ˇı ́kov ́a 等，2007；Alkema 等，2007；Radtke 等，2002；Fuentes 和 Raftery，2005）。

6 其他方法和目标

在这里，我们简要概述了其他代理建模和下游任务。

6.1 回归树

在某些情况下，模拟器均值 M 可能具有不连续性或 “状态变化” ，其中输入空间的一部分与另一部分相比，y 和 x 之间存在非常不同的关系（即非平稳性）。回归树 (Breiman 等, 1984) 形成了一类在这些情况下很有用的方法。它们在某些输入是分类而不是数字的情况下也很有用。通过将输入空间划分为相互排斥的区域来处理这些问题，在这些区域中适合独立的代理（GP 或其他回归方法）。

两种方法：树状 GP (TGP Gramacy and Lee, 2008) 和贝叶斯加性回归树 (BART Chipman 等, 2010) 已得到广泛应用。两者都使用数据自动划分输入空间，依赖贝叶斯计算，都有公共软件：TGP in tgp on CRAN (Gramacy and Taddy, 2016; Gramacy, 2007)；几个 R 包中的 BART，包括 BayesTree（Chipman 和 McCulloch，2016）和 BART（McCulloch 等，2019）。

Rulli`ere 等的其他方法。 (2018)，以及通过 Voronoi 镶嵌而不是树木（例如，Kim 等，2005；Rushdi 等，2017；Park 和 Apley，2018），受到的关注较少。普拉托拉等。 (2020) 通过将 M 建模为贝叶斯回归树的总和（如 BART）和作为贝叶斯回归树的乘积的内在方差 \sigma^2_v(x)，采用类似于在 第 3.2 节 中。

Konomi 等探讨了利用 KOH 方法和使用 TGP（ 第 6.1 节 ）中的回归树的校准方法。 (2017)。在分区的每个终端节点中，为计算机模型输出假设一个具有独立常数固有方差项的 GP。还为差异项部署了独立的 GP。尽管 \sigma^2_v 在每个终端节点都是常数，但常数在终端节点之间可能会有所不同，因此会自动包含异方差性。

6.2 定性输入

分类（定性）变量通常出现在随机模拟器中，尤其是那些包含人类行为特征的模拟器。虽然回归树能够处理分类输入（Broderick 和 Gramacy，2011；Gramacy 和 Taddy，2010），但 GP 作为平滑模拟器输出的替代品可能更有效。

钱等。 (2008)，周等 (2011), 和 Chen 等 (2013) 描述了扩展用于数值输入的内核以合并定性变量的方法。粗略地画画，他们的方法将两个输出 y(x_i) 和 y(x_j) 之间的相关性作为两个相关函数的乘积：处理连续输入 w 的 Cc(wi, w_j) 和 Cq(zi, zj) 对于定性变量 z。构建 Cq 的一种简单方法采用 Cq(wi, w_j) = K \infty k=1 τk,wik,w_jk \tag{6.1}

其中 K 是定性变量的数量，τk,wik,w_jk 表示 wik 和 w_jk 之间的相关性。 τk,wik,w_jk 的一个例子是： τj,wik,w_jk = \exp{-(φik + φjk)I[wik 6= w_jk]} \tag{6.2}

其中 I 是指示函数（如果参数为真则 = 1，如果为假则 = 0），并且 φ ^{\top} 0。引用的参考文献还提供了其他建模 τk,wik,w_jk 的方法。存在替代方法，例如 Zhang 等 (2018) 将潜在变量用于定性模型。

6.3 优化

一个常见的实验目标是最大化模拟器的输出，即找到最大化输出 y(x) 的输入 xmax。为了最小化，将 y(x) 替换为 -y(x)。优化通常是一个顺序过程，其中选择连续的 xs 以越来越接近最优 xmax——一个顺序设计问题（参见 第 4 节 ）。对于随机模拟器，y(x) 是随机的，最优值的定义不太具体——每次模拟器在相同的 x 下运行时输出都是不同的。因此，兴趣通常在于最大化感兴趣的非随机数量，例如平均值 M 或可能是另一个标量，例如第 q 分位数。

对于确定性模拟器，贝叶斯优化（Mockus 等，1978；Jones 等，1998）是一种流行的技术。一组初始运行用于构建 GP 代理，并通过最大化 “获取函数” α(x) 选择新运行。迭代选择 x_{new} = arg maxx α(x)) 提供了一个逐渐改进的最大值估计。 α(x) 的一个广泛使用的选择是预期改进 (EI)：αEI(x) = E[max (y(x) - ymax, 0)]。 \tag{6.3}

最大化 EI 选择输入 x_{new} 最大化已观测运行的最大值 ymax 的预期增长。 y 由 GP 建模： ) \tag{6.4}

其中 \mu_N (x) 是 GP 的预测平均值，σN (x) 是其标准差，φ 是标准正态密度，Φ 是标准正态分布函数。

近年来，替代获取函数在贝叶斯优化方面产生了广泛的工作，主要是在机器学习文献中。改进概率 (Kushner, 1964) 是一个早期示例，其他示例，例如 GP 置信上限 (GP-UCB)（Srinivas 等，2009），考虑同方差模拟器误差。最近的摘要可以在 Frazier (2018) 中找到。”

Picheny 等讨论并比较了具有恒定固有噪声的随机问题的替代标准。 (2013);用上述方法称为 “插件” 方法。 DiceOptim（Picheny 等（2016）；Picheny 和 Ginsbourger（2014））中提供了一个用于实现其中几个选择的 R 包。贾拉利等 (2017) 也对异方差噪声进行了类似的比较。

水平集估计的相关目标是找到输出超过阈值 T 的区域，也可以使用类似于 EI 的顺序标准作为目标。一个简单的标准是最大轮廓不确定性 (MCU)，其中新点是根据它们被认为与 T 的接近程度以及该点的不确定性程度的加权和来选择的。吕等。 (2018) 在这里提供一些讨论。 hetGP 中也实现了这种方法。

使用高斯过程进行优化，特别是在存在潜在异方差（并且可能是非正态）的内在可变性的情况下，是一个有趣的研究问题，可能值得对其进行调研。尽管如此，此处提供的参考资料应该提供了很好的介绍。

6.4 敏感性分析

确定和测量输入对输出的影响通常是任何模拟器实验的一部分。这样做有助于科学地理解系统，并能够筛选出潜在的多余变量。这个目标有很多相关的名称：敏感性分析、筛选、变量选择等，但总体目标大体相同——总结和衡量每个输入的影响。

对于确定性模拟器，Sobol 指数（Sobol，1993）被广泛使用。在模拟器的输入上假定概率分布，以表示它们的变化范围。然后，函数方差分析 (ANOVA) 分解将模拟器输出的变化分成多个分量，每个分量代表输入变量 x_j 或输入变量组合的单独贡献。然后将 Sobol 指数计算为由组件解释的总模拟器输出变化的百分比。关键 Sobol 指数包括主效应（单独由单个 x_js 解释的变异百分比）和由与其他输入的交互加性效应解释的变异。计算组件需要大量运行，但使用代理 GP 使计算可行（Schonlau 和 Welch，2006；Marrel 等，2009）。 Oakley 和 O’Hagan (2004) 对敏感性进行了全面的讨论

Marrel 等的两个扩展。 (2012) 和 Hart 等 (2017) 的 Sobol indices for stochastic simulators 为随机模拟器产生以下表达式： y(x) = y(x, seed) \tag{6.5}

其中 x 是可控输入集。输入种子负责输出随机性，代表内在可变性，有时称为种子变量。与确定性模拟器一样，假设概率分布（通常是均匀的）代表可控输入的变化范围。

在 Marrel 等 (2012)，通过功能方差分析分解分析随机模拟器均值的总变化，并根据每个组件解释的总模拟器变化的百分比计算 Sobol 指数。也可以计算由种子变量 seed 解释的变异，表示由内在方差解释的总变异。此外，可以单独进行内在方差 \sigma^2_N (x) 的敏感性分析，以收集有关哪些输入变量对异方差影响最大的信息。

哈特等。 (2017) 假设模拟器可以在具有相同种子的不同输入 x 下运行。如 第 3.2 节 所述，他们没有为均值和方差构建一个联合随机模拟器代理，而是为多个种子构建一个单独的代理。对于每个种子，他们获得每个 Sobol 指数的实现，并通过汇总实现，他们获得指数的分布。

关于模型选择的广泛文献可能有对随机模拟器有效的对应物。但是即使对于确定性模拟器来说，完全令人满意的方法仍然有些难以捉摸。

7 结束语

从这篇综述中可以得出几个关键信息，每个信息都指向开放或新的研究问题：

高斯过程代理。 GP 之所以被广泛讨论，是因为它们提供了一种灵活的方式，允许数据告知底层流程的形状。此外，它们可以是不确定性的有效预测因子和量化因子。诊断随机模拟器（在确定性设置中可用（Bastos 和 O’Hagan，2009））的 GP 中的缺点尚未确定。

如 第 3 节 所述，神经网络（深度学习）方法正在积极使用和研究中，其中一些方法可能与 GP 相结合，提供有前途的研究方向（Schultz 和 Sokolov，2018）。

此外，可能难以有效地捕获非正态变异性。用尽可能少的模拟来做到这一点，同时也适当地量化各种不确定性，可能是随机模拟器分析的一个重要研究方向。更广泛的分位数回归文献可能是一个很好的起点。

设计。随机模拟器不同于确定性模拟器，因为它们需要更大的样本量并允许使用重复，其处理通常是临时的。这引出了 第 4 节 提出的问题，形成了一个重要的研究方向。设计尺寸的经验法则即使不完美也很有用，适用于确定性模拟器 (Loeppky 等, 2009b)，但缺乏随机模拟器。

校准。考虑校准中的模型差异至关重要，但没有明显的 “一刀切” 方法。需要进行广泛的实证比较，并指导哪些策略在哪些条件下有效。评估不同方法的有效性可能具有挑战性（参见 McKinley 等，2018，关于 ABC 和 HM 之间的一项比较），但这是迫切需要的。

模拟器复杂性对于复杂的随机模拟器，获得足够的运行可能是不可行的。在某些情况下，模拟器可以替换为一个不太复杂的模拟器（例如 Molina 等，2005），它可以捕获关键特征并允许进行足够数量的模拟。已经探索了另一种途径，将随机模拟器与确定性模拟器耦合（Baker 等，2020）作为处理低模拟预算的方法。多保真建模，其中多个不同复杂性的模拟器耦合在一起（Kennedy 和 O’Hagan，2000；Kennedy 等，2020）是一个有前途的解决方案

同样，某些输出可能比其他输出噪声小，噪声较小输出的建模可以改进噪声较大输出的建模。例如 Wang 和 Ng (2020) 使用模拟器的期望来改进噪声分位数的估计。这与历史悠久的更广泛的方差减少文献有关 (Barton 等, 2017)。方差减少已应用于许多示例中，但它在 ABM 中的使用并不明显，这可能是由于 ABM 中随机元素的大量使用。固定随机模拟器中的初始种子在敏感性分析中发挥了作用（参见 第 6.4 节 ），但将有关内在随机性的信息用于更广泛的目的是一个悬而未决的问题。

这篇综述致力于提高人们对处理随机模拟器的现有工具和策略的认识，并为有兴趣利用最新统计方法的从业者提供一个起点。尽管问题普遍存在且具有挑战性，但该领域的统计研究仍然不足。这些问题提出了计算和技术问题，以及理论和哲学问题。目前的解决方案往往是有能力的，但缺乏对不同解决方案之间的全面比较，也缺乏对它们对复杂情况（例如非常大的数据集）的普遍性的测试。希望该评论能激发统计研究人员参与讨论的开放性问题。