论文笔记 计算机模型的贝叶斯校准 Bayesian calibration of computer models

计算机模型的贝叶斯校准 Bayesian calibration of computer models

摘要

我们考虑对使用复杂数学模型进行近似的系统进行预测和不确定性分析。此类模型作为计算机代码实现，通常是通用的，因为通过适当选择模型的一些输入参数，代码可用于预测系统在各种特定应用中的行为。.然而，在任何特定应用中，必要参数的值可能是未知的。在这种情况下，系统在特定环境中的物理观察用于了解未知参数。通过调整参数将模型拟合到观测数据的过程称为校准。校准通常通过特设拟合来实现，并且在校准之后，使用带有拟合输入值的模型来预测系统的未来行为。我们提出了一种贝叶斯校准技术，它在两个方面改进了这种传统方法。首先，预测允许所有不确定性来源，包括拟合参数的剩余不确定性。其次，他们试图纠正模型的任何不足之处，这些不足之处是观察到的数据与模型预测之间的差异所揭示的,即使是最佳拟合参数值也是如此。该方法通过使用来自托木斯克核辐射释放和更复杂的模拟核事故演习的数据来说明。

关键词：校准；计算机实验；确定性模型；高斯过程；插值；模型不足；敏感性分析;不确定性分析

1.概述

1.1 计算机模型和校准

各种科学使用数学模型来描述原本很难分析的过程，这些模型通常在计算机代码中实现。 通常，数学模型非常复杂，生成的计算机代码很大，并且就单次运行所需的计算机时间而言可能很昂贵。 然而，运行计算机模型将比直接观察过程便宜得多。 Sacks、Welch、Mitchell 和 Wynn (1989) 给出了几个例子。 我们认为的代码是确定性的，即:使用相同的输入运行代码总是会产生相同的输出。

计算机模型通常被设计为适用于广泛的特定环境。 但是，要使用模型在特定环境中进行预测，可能需要首先使用一些观察到的数据来校准模型。 为了说明这个过程，我们介绍一个简单的例子。 2.2 节详细描述了另外两个示例。

**********************(药物的性质决定某些转移比率，从而决定探索后果)********************** _____________校准是调整未知速率参数直到模型的输出适合观察到的数据的活动___________

为了决定新药的剂量方案（例如片剂的大小、频率和释放速率），使用了药代动力学模型。这模拟了药物通过患者身体的各个“隔室”的运动及其最终消除（例如通过化学反应或排泄）。这样的模型允许探索任何给定剂量方案的后果。然而，要将模型用于特定药物，必须指定药物在不同身体部位之间移动的速率，例如从胃转移到血液的速率，或从肝脏中消除的速率。这些比率中的一些或全部将特定于所讨论的药物。正是通过要求用户将它们指定为代码的输入，该模型才实现了对各种药物的适用性。当然，对于给定的药物，这些比率是未知的，因此进行实验以获得观察数据。无法自行获取费率。相反，数据是对药代动力学模型的某些输出的观察，例如在某些时间点血液或尿液中的浓度。粗略地说，校准是调整未知速率参数直到模型的输出适合观察到的数据的活动。

更一般地，计算机模型将具有许多特定于上下文的输入，这些输入定义了要使用模型的特定情况。通常情况下，当一个或多个特定于上下文的输入的值未知时，可以使用观察来了解它们，这是校准。

在当前实践中，校准总是包括搜索一组未知输入的值，以使观察到的数据在某种意义上尽可能接近模型的相应输出。这些值被视为特定于上下文的输入的估计值，然后通过将这些输入设置为它们的估计值，该模型用于预测该上下文中过程的行为。

显然，这种“插件”预测将特定于上下文的输入视为已知的。现实情况是它们只是被估计的，关于这些输入的剩余不确定性应该在模型的后续预测中得到承认。

我们在本文中提出了一种贝叶斯方法来校准计算机模型。我们将未知输入表示为参数向量 θ 。使用观察到的数据，我们得出θ 的后验分布，它特别量化了关于 θ 的“剩余不确定性”。随后使用计算机模型进行预测时，这种不确定性得到了充分考虑，通过获得预测分布。 Aitchison 和 Dunsmore (1975) 阐述了贝叶斯预测推理的原理，但计算机代码校准问题具有一些复杂的特征。

我们的方法将计算机代码视为简单的“黑匣子”。我们不使用有关代码实现的数学模型的信息，除非这可以通过有关输入和输出之间关系的先验信息来表示。显然，打开黑匣子并利用其结构的方法可能被证明比我们的方法更强大，但它们的应用也会相应地更复杂。这是未来研究的一个潜在重要课题。

1.2 本文大纲

第 2 章详细分析了与计算机代码相关的问题。 特别是，除了上下文特定输入的不确定性之外，计算机代码的使用还有其他几个不确定性来源，而用于预测的校准只是计算机代码输出统计分析中的几个主题之一。 第 2 章首先仔细研究了计算机模型中的不确定性，并通过一些详细示例说明了可能的来源。 最后回顾了以前关于校准和相关问题的工作。

我们的贝叶斯方法基于通用贝叶斯方法来推断未知函数。 第 3 章回顾了相关理论和已发表的工作。第 4 章介绍了我们的基本模型和分析，第 5 章讨论了实际应用问题。第 6 章介绍了基于真实数据的案例研究。 第 7 章提供了一些结论和进一步工作的方向。

2.计算机代码输出的统计分析

2.1 计算机模型的不确定性

计算机模型的广泛应用伴随着对量化其使用中普遍存在的不确定性的广泛关注。 以下是对各种不确定性来源进行分类的一种方法。

不确定性	来源
模型不足	没有模型是完美的
参数不确定性	计算机代码输入值视为模型的未知参数
残余变异性	由于无法识别的条件导致
参数可变性	输入中指定的某些条件不受控制和未指定
观察误差	观察错误
代码不确定性	给定任何特定输入配置的计算机代码的输出实际上是未知的

2.2 例子

能够将本文中的各种概念、定义和符号与实际计算机代码的一些具体示例联系起来会很有帮助。以下只是可以给出的众多示例中的两个。

2.2.1 高斯羽流模型

在放射防护领域，一个简单的高斯羽流模型（Clarke，1979 年）用于预测放射性物质在意外释放后的扩散和随后的沉积。在这些情况下，运行更复杂模型所需的详细输入信息不可用。

代码输入可分为定义事故发生时的大气条件（风向、风速和大气稳定性）和定义释放性质的输入（源项、源位置、释放高度、释放持续时间和沉积）速度。（条件决定性质输入）

放射性核素的扩散是一个高度复杂的过程，涉及各种无法直接观察到的化学和环境过程。在高斯羽流模型中进行了许多简化假设，通常会导致模型高度不充分。例如，假设风的速度和方向在释放的粒子的传播时间内保持不变。

即使对于这个简化的模型，许多输入也具有与之相关的参数不确定性。源项表示释放的物质量，沉积速度表示地面空气中物质沉积在地面上的速率，这些都是示例。沉积速度的合适值很难确定（参见 Jones (1981)）。默认值通常基于粒子的大小和材质通过的地形类型来使用。释放的高度和持续时间，以及风速和风向也可能具有相关的参数不确定性。高斯羽流模型很便宜。我们可以在很短的时间内进行数千次运行，因此出于实际目的，代码不确定性不是问题，我们可以处理已知的功能。这不适用于更复杂的大气扩散模型。
在实际发布过程中，不同组织使用各种设备进行沉积测量。进行这些测量的过程并不简单。通常从现场采集草样并在实验室进行分析。在测量过程的每个阶段都会引入测量误差。

该讨论涉及在特定事故场景中使用羽流模型，其中参数不确定性表示对实际释放的输入参数真实值的置信度。我们可能还想考虑风险评估环境，我们想预测未来可能的污染。事故随时可能发生，因此源项、风速和风向是随机的，输入受参数可变性的影响。

2.2.2 水文模型

水文模型用于预测地下水流量，例如预测土壤中污染物的移动或降雨后河流的排放量。我们考虑 Romanowicz 等人描述的模型。 （1994 年），它预测流放电。模型的输入包括测量的降雨数据的时间序列、蒸散速率、剖面刚饱和时土壤的平均有效透过率 TO 以及与地下排水水平相关的常数 m。

关于 TO 和 m 的参数不确定性源于该模型用于预测水流通过无法直接观察到的复杂地质结构的事实。测量的降雨数据将包括随机测量误差，这是参数可变性的一个例子。模型的不足源于建模过程引入的简化，如高斯羽流模型。

Romanowicz 等人。 (1994) 使用实际流量测量来了解 m 和 TO 的值。然而，正如他们所指出的，不可能测量流动过程的真实价值，因为我们只能对流动的异构模式进行点测量。这与第 2.1 节中描述的剩余可变性的概念有关，我们将“真实流量”定义为平均值。

2.3. 统计方法和以前的工作

计算机代码输出统计分析领域的一些早期工作主要涉及插值，即给定包含输入配置样本输出的数据，问题是估计代码在其他输入配置下的输出尚未运行。当代码特别大且运行成本特别高时，这是相关的。 Sacks、Welch、Mitchell 和 Wynn (1989) 对这项工作进行了重要回顾，最近的一些参考文献是 Currin 等人。 （1991），莫里斯等人。 （1993 年），贝茨等人。 (1995) 和肯尼迪和奥哈根 (2OOOa)。

在这项工作中考虑的唯一形式的不确定性是代码不确定性。模型不充分、残余变异和观察误差不相关，因为没有尝试预测与计算机模型输出不同的真实过程，并且没有使用对该过程的观察。所有输入参数都假定是已知的，并且不受参数不确定性或参数变化的影响。这项工作中使用的统计方法基于将计算机代码输出表示为其输入的未知函数，并将该函数建模为随机过程。

统计方法已经解决了相当长一段时间的另一个问题是不确定性分析。不确定性分析的目的是研究由输入概率分布引起的代码输出分布。输入参数分布可以是参数不确定性的公式，即其值未知的参数，或参数可变性的公式，即其值未指定的参数。最简单的不确定性分析方法是蒙特卡罗解决方案，其中输入的配置是从它们的分布中随机抽取的。然后针对每个样本输入配置运行代码，生成的输出集是来自要研究的输出分布的随机样本。参见 Helton (1993) 的评论。

用于不确定性分析的蒙特卡罗方法很简单，但当代码运行成本很高时变得不切实际，因为需要大量运行。与简单的蒙特卡罗抽样相比，拉丁超立方抽样（McKay 等人，1979；Stein，1987；Owen，1992）声称效率更高；参见例如 Crick 等人。 （1988 年）和 Helton 等人。 (1991)。阿斯莱特等人。 (1998) 将蒙特卡洛方法与代码的统计插值结合起来进行不确定性分析，有效地使用插值器作为代码的廉价替代品。

Haylock 和 O’Hagan (1996) 提出了一种完全不同的基于高斯过程先验模型的不确定性分析贝叶斯方法。他们推导出输出分布的后验均值和方差，并将其推广到 Oakley 和 O’Hagan (1998) 对输出分布的分布函数和密度函数的后验估计。

除了代码不确定性之外，这些不确定性分析方法还考虑了参数不确定性和参数变化。但是，目标仍然集中在代码输出上，在这种情况下是以输出分布的形式，而不是过程本身。因此，没有对模型不足、残差或观测误差进行处理。

与计算机代码相关的另一个问题是敏感性分析，其目标是表征代码输出如何响应输入的变化，特别是涉及识别输出相对敏感或不敏感的输入。有关该主题的大量文献的一个很好的来源是 Saltelli 等人。 (2OOO)。尽管其中大部分没有使用统计方法，但也有一些值得注意的例外。参见例如 Helton (1993)、Morris (1991)、Iman 和 Conover (198O)、Welch 等人。 (1992)、Morris (1991)、Homma 和 Saltelli (1996) 以及 O’Hagan 等人。 (1999)。在德雷珀等人。 （1999 年），敏感性分析被应用于一系列“场景”，使用 Draper（1995 年）的模型平均思想。与插值一样，这些敏感性分析的统计方法仅考虑计算机代码输出的所有统计分析所共有的不确定性来源，即代码不确定性。

我们对该领域调查的最后一个主题是校准。如 1.1 节所述，估计未知参数的传统方法是通过临时搜索最佳拟合值。因此，一些考虑因素是观察误差、残差和模型不足，但只是通过测量拟合差异来隐含地进行。该度量通常不是通过以任何显式方式对这些误差项建模而开发的，并且通常完全是启发式的。此外，由于估计值随后被视为已知，因此后续预测不考虑（剩余）参数不确定性。

相比之下， Romanowicz 等人的广义似然不确定性估计方法。 (1994) 确实充分考虑了参数的不确定性。该方法实际上是贝叶斯方法。从未知输入的先验分布中抽取初始蒙特卡洛样本，然后通过似然项对其进行加权。使用所有采样的输入配置进行预测，结果是来自后验预测分布的加权样本。例如，加权样本均值是对该预测分布均值的估计。

在广义似然不确定性估计方法中，可以允许仅具有运行样本的代码不确定性，并且还可以通过扩展蒙特卡罗方法来允许参数变化以允许在预测时从未指定的输入中抽取阶段，但这些在迄今为止发表的文献中都没有完成。与更传统的校准方法一样，可能性是相当启发式的，并且只是非常松散地基于对代码输出和实际过程之间的差异进行建模。模型不充分、残差和观测误差没有被明确区分或建模。在预测中没有考虑它们，目标仍然是估计代码输出而不是现实。

Craig 等人给出了另一种贝叶斯校准方法。 (1996, 2OO1)。他们对代码输入和输出之间的关系进行了建模，类似于在该领域的许多其他工作中提到的高斯过程模型。在克雷格等人。 (1996)，主要目标是使寻找最佳拟合校准更加有效和系统化。他们的建模反映了搜索的迭代性质，他们的方法采用了 Goldstein (1986, 1988) 和 Wooff (1992) 的贝叶斯线性哲学，而不是完全指定的贝叶斯分析。该方法在 Craig 等人中得到了扩展。 (2OO1) 处理校准后的预测问题。

考克斯等人。 (1992) 描述了一种校准方法，该方法类似于传统的最佳拟合参数搜索，但它用使用高斯过程模型获得的更便宜的插值器代替了昂贵的代码。他们的方法没有考虑预测阶段剩余的参数不确定性。另见 Cox 等人。 (1996)。

Raftery 等人尝试结合先前专家对校准参数和模型输出的意见。 （1995）使用他们称为贝叶斯综合的方法。这受到 Wolpert (1995) 和 Schweder 和 Hjort (1996) 的批评，并在后续论文中 Poole 和 Raftery (1998) 提出了一种称为贝叶斯融合的替代方案。这两种方法都没有明确地认识到模型的不足，并且底层的计算机代码应该足够简单，可以忽略代码的不确定性。

我们在本文中的分析是第一次尝试对计算机模型的校准和后续使用中出现的所有不确定性来源进行建模并明确考虑。我们坦率地承认，我们目前实施的方法不是完全贝叶斯方法，因为我们通过（近似）后验模式估计超参数。充分考虑超参数中的不确定性是进一步研究的另一个主题，但我们相信我们目前的方法是迄今为止对计算机代码不确定性最全面的处理，我们认为它在实践中可能是足够的。

3. 函数的贝叶斯推理

3.1 高斯过程

高斯过程的使用在2.3节被提高过几次，令 $f(.)$ 为输入 $x$ 映射到 $y=f(x)$ 的映射关系，其中 $x\in \chi$ ，$\chi$可以是任意的，以最典型的 $R^q$ 为例，$x=(x_1,x_2,...x_q)$，我们将 $y$ 视为标量，$y\in R$，也可以视为向量$y\in R^{q^{\prime }}$ 。

从形式上看， $f(.)$ 可表示 $f(X_1),f(X_2),f(X_n)$ ,且均为正态分布。
这个分布的均值为 $m(.)$，$m(x)=E[f(x)]$，协方差 $c(.,.)$,$c(x_1,x_2)=cov[f(x_1),f(x_2)]$,那么
$f(.)$~$N${ $m(.),c(.,.)$ }

高斯分布是一种灵活且方便的分布函数。

3.2. 建模问题

在实践中使用高斯过程的原因与在统计理论和建模中普遍使用正态分布的原因大致相同。 它们方便、灵活并且通常非常现实。 当然，重要的是正态性，特别是联合正态性，是对 $f(.)$ 的先验知识或信念的合理表示。 变换在这种情况下可能很有用，就像它们在正常理论模型的其他应用中一样。

鉴于高斯过程是一个合理的建模选择，然后应指定均值和协方差函数以反映有关 $f(.)$ 的详细先验知识。 例如，如果平稳性是先验置信度的一个特征，那么对于任何 $d\in \chi$，$f(x)$ 的先验分布与 $f(x+d)$ 的先验分布相同（并且定义了对 $\chi$ 的加法运算） , $m(.)$ 将是一个常数，而 $c(x,x^{\prime })$ 只是 $x-x^{\prime }$ 的函数。

通常，$m(.)$ 可以是 $\chi$ 上的任何函数，但 $c(.,.)$ 必须是对于每个$n=1,2,...$ 对应的 $f(x_1),f(x_2).....,f(x_n)$, 的协方差矩阵。包括元素 $c(x_i,x_j)$ 对于所有 $x_1,...x_n\in \chi$ 。 都是非负定的。 Cressie (1991) 给出了协方差函数在这个意义上为非负的一些条件。

一个好用的方法是对 $m$ 和 $c(.,.)$ 分层 建模。 对于 $m(.)$ 。 我们可以使用线性模型结构

$m(.)=h(.{)}^T\beta$ （1）

其中 $h(.)=(h_1(.),h_2(.)...h_p(.){)}^T$（已知）， $\beta (.)=({\beta }_1,{\beta }_2...{\beta }_p{)}^T$（未知），因此 $h(.)$ 描述了一类形状并且模型 (1) 可以用这个类中的一个函数来近似 f $(.)$ 。比如，当x是标量时， $h(x)=(1,x...x^{p-1}{)}^T$, 表示 $m(.)$

作为 $\beta$ 的先验分布，多元正态分布是一个方便的选择。例如，它具有 $f(.)$在积分出 $\beta$ 后仍为高斯过程的性质。当然，应指定先验分布以反映真实的分布而不是方便，但在实践中，关于超参数例如 β的先验信息通常会很弱。因此，经常使用通过不适当的均匀密度$p(\beta )\propto 1$来表示弱先验信息的常规表示。
请注意，我们可以通过编写来分离均值和协方差结构

$f(x)=m(x)+e(x)=h(.{)}^T\beta +e(x)$ （2）

使用模型 (1)，其中 $e(x)$ 是零均值高斯过程，具有协方差函数$c(.,.)$。 $c(.,.)$ 的建模非常重要，因为通过$f(x)$ 和 $f(x_t)$ 之间的相关性，我们表达了一个观点，即如果 $x$ 和 $x_t$ 在 $\chi$中足够接近那么 $f(x)$ 和 $f(x_t)$ 应该相似 ，因此表达了对 $f(.)$ 平滑度。 在接下来的部分中，我们一般将采用具有第一阶段的 $c(.,.)$ 的层次模型

$c(.,.)={\sigma }^{2}r(x-x^{\prime })$ （3）

其中 $r(.)$ 是具有属性 $r(0)=1$ 的相关函数。这个公式表达了关于 $e(.)$ 的先验信息的平稳性； 我们有一个共同的（未知的）方差 ${\sigma }^2$ 和仅取决于$x-x^{\prime }$ 的相关性。 相关函数进一步用其他未知超参数表示； 例如，如果 $x=(x_1,...,x_q)$ 是一个向量，

$r(x-x^{\prime })=exp(-{\sum }_{j=1}^q{w}_j(x_j-x_j^{\prime }{)}^2)$ (4)
然后在下一阶段，我们将表示 方差 ${\sigma }^2$和 粗糙度参数 $w_1，.....,w_q$. 的先验分布.

当然等式（4）只是一种可能的公式。 一个已被广泛使用的更通用的表达式将 $(x-x_t{)}^2$替换为 $||x-x_t|{|}^a$，其中 a 可能具有指定值或另一个超参数。 更一般地说，我们可以在每个维度上允许不同的 a 。 等式 (4) 的另一个推广是设置$r(d)=exp(-d^T\Omega d)$，其中Ω是未知对称正定矩阵，在等式 (4) 中具有 $\Omega =diag（w_1,\dots w_q)$ 的形式。
在地质统计学中，通常会花费大量精力来估计 ${\sigma }^{2}r(x-x^{\prime })$ ，或者等价于 σ 2 ( r ( 0 ) − r ( x − x ′ ) ) σ^2 (r(0)-r (x- x'))