【摘 要】 当标准似然估计方法在计算上不可行时，我们建议使用深度学习来估计统计模型中的参数。我们展示了如何从最大稳定过程中估计参数，其中即使使用小数据集进行推理也非常具有挑战性，但模拟很简单。我们使用来自模型模拟的数据作为输入并训练深度神经网络来学习统计参数。我们基于神经网络的方法为当前方法提供了一种有竞争力的替代方法，这一点在相当大的准确性和计算时间改进中得到了证明。它作为统计参数估计中深度学习的概念证明，可以扩展到其他估计问题。

【原 文】 Lenzi, A. et al. (2021) ‘Neural Networks for Parameter Estimation in Intractable Models’. arXiv. Available at: http://arxiv.org/abs/2107.14346 (Accessed: 15 November 2022).

1 简介

由于数据存储和传感器技术的进步，越来越多的大型数据集不可避免地显示出复杂的依赖关系，给统计建模和预测带来了新的机遇和新的挑战。环境过程给出了一个重要的例子，数据集可以在多个空间和时间变化尺度上具有复杂的相互作用。具体而言，为考虑时空依赖性的非高斯数据设计的参数统计模型的推断在计算上具有挑战性。例如，考虑多变量极端模型，其中数据通常在大量空间位置进行采样，但计算复杂性大限制了它们的实用性。为了规避这个瓶颈，已经开发了避免全似然计算的模型近似。一种这样的近似基于复合似然，通常由成对的观察结果构成（Tawn，1988 年；Padoan 等人，2010 年；Ribatet 等人，2012 年）。然而，即使最大成对似然估计量在温和的正则性条件下是强一致且渐近正态的（Xu 和 Reid，2011），但不存在用于最优选择复合似然项的明显策略，效率的损失使得高维推理成为可能数据仍然具有挑战性（Varin 等人，2011 年；Huser 等人，2016 年）。

作为空间极值模型中复合似然推理的替代方法，Stephenson 和 Tawn (2005) 建议用它们的出现时间来增加分量块最大值数据，并表明根据最大值是否重合来划分位置可以大大简化似然函数。为了避免修复这些分区时出现的偏差，Huser 等人。 (2019) 通过将分区视为潜在变量，构建了一个随机期望最大化算法，用于最大稳定过程的精确最大似然估计。同样，在贝叶斯设置中，Thibaud 等人。 (2016) 开发了马尔可夫链蒙特卡罗算法，用于以分区值为条件的最大稳定过程。这些方法的一个可能问题是它们假设时间独立，并且这种时间依赖性可能导致多少偏差尚不清楚。此外，基于似然的技术仍然受到评估多元最大稳定累积分布函数的计算瓶颈的阻碍。依赖模型的一个关键方面是模拟通常是快速且易于处理的，尽管似然计算很繁琐。利用这一点，Erhardt 和 Smith（2012 年）采用了一种基于模拟的方法，使用近似贝叶斯计算 (ABC) 来估计最大稳定模型的参数。在其最基本的形式中，ABC 首先从先验分布中采样一组参数，然后将其用于在假设的统计模型下模拟数据。对于每个新数据集，必须重复评估模拟和观察数据之间差异的过程。因此，前面提到的方法都不适合大型数据集。因此，为计算上具有挑战性的可能性开发新的推理方法仍然是一个开放的挑战，就像最大稳定过程和许多其他应用程序的情况一样。

在这项工作中，我们提出了一种使用深度神经网络 (NN) 来估计棘手统计模型参数的方法。深度学习算法利用神经网络作为多个互连层来预测或分类数据集之间的复杂映射，近年来它们在预测高维非线性过程方面取得了显著成功。关于直接应用于参数估计的深度神经网络技术的文献很少。在 ABC 框架内，Creel (2017) 提议使用信息统计来训练深度神经网络，经典或贝叶斯间接推理可能以此为基础。姜等。 (2017) 根据模拟数据预测汇总统计数据，然后将其用作 ABC 程序后验均值的估计。最近，Radev 等人。 (2020) 提议使用来自数据的信息汇总统计作为可逆神经网络的输入来执行贝叶斯推理。然而，这些论文没有提供关于如何在实践中为网络构建合适的汇总统计的一般指导。通过使用传感器测量（Liu 等人，2018 年）、时间序列数据（Cremanns 和 Roos，2017 年）和空间数据（Gerber 和 Nychka， 2020 年）。鲁迪等。 (2021) 提出使用 NN 中的密集层和卷积层来解决确定性逆问题，其中来自常微分方程的非线性系统的参数是根据模拟数据估计的。据我们所知，目前的工作是第一项利用深度神经网络方法的研究，目的是仅基于模拟数据和参数进行推理。

与 ABC 类似，我们基于模拟执行统计推断，而不是直接评估计算量大或难处理的似然函数。与 ABC 的一个关键区别在于，我们的方法不是将相似性作为模拟数据与实际观察的信息汇总统计的函数来衡量，而是学习直接在数据和参数之间执行映射。因此，我们的方法可以看作是自动 ABC 的一种形式，其中数据压缩是在神经网络内部完成的，而没有选择潜在的高度敏感指标的困难。我们的新结构使用来自统计模型的模拟数据作为深度神经网络的输入。输出由用于模拟输入数据的参数组成。我们利用可行性来模拟任意大量的训练数据，以确保深度神经网络尽可能接近实际参数。深度神经网络训练完成后，将观测数据作为测试数据的输入，返回感兴趣的模型参数作为输出。我们表明，与传统（近似）似然方法相比，深度神经网络可以以更高的精度估计空间极值的最大稳定模型的参数，并显着加快计算速度。随着参数数量的增加，ABC 方法的计算效率会大幅下降，我们提出的框架也可能会出现同样的情况。这项工作在参数估计方面取得了可喜的成果，并作为在传统技术失败时机器学习如何用于统计推断的概念证明。

我们在从最大稳定模型确定参数的问题上展示了我们的框架。极端的空间数据可以被视为 2D 图像的特例，这使得具有卷积层的神经网络成为有前途的候选者，因为它们在学习复杂图像的特征方面取得了成功。在这项工作中，我们首先对 Brown-Resnick（Kabluchko 等人，2009 年）和 Schalther 的（Schlather，2002 年）模型进行了模拟示例，并将我们的方法的准确性与标准成对似然技术进行了比较。我们表明，CNN 架构可以快速适应，并且可以提供比典型方法提供的更少偏差和更少可变的预测。然后，使用美国中西部 10 天的再分析温度数据块，我们说明了所提出的框架在准确性和计算节省方面优于标准方法的性能和优势。我们的方法相对于成对似然估计方法的另一个优势是不需要该领域的详细知识（例如，极端），因为目标统计计算是在深度神经网络内部完成的。因此，我们基于深度神经网络的框架不需要对数据生成过程进行任何假设。相反，它被设计成具有灵活且计算上吸引人的形式，通过训练学习网络的权重，从而可以准确表示一大类函数。虽然我们的重点是空间极值建模，但只要可以从概率模型进行模拟，新方法就可以用于其他推理问题。

本文的其余部分安排如下。第 2 节概述了所提出的方法并简要介绍了 CNN。在第 3 节中，我们回顾了基于成对似然的最大稳定过程和推理的定义，并进行了模拟研究以评估估计器的性能。在第 4 节中，我们在温度最小值数据集上说明了我们的方法，并在第 5 节中进行了讨论。

2 方法论

本节概述了我们使用深度神经网络估计统计模型参数的新框架，并简要概述了示例中使用的网络类型。

2.1 参数估计框架

在下文中，我们将使用术语参数来指代感兴趣的统计模型中的参数；术语 NN 权重表示 NN 的参数。与经典统计推断一样，我们的目标是使用观察到的数据通过估计假设概率模型的参数来推断随机过程的潜在特征。假设 Y 是一个随机变量，其值位于可测量空间中，概率密度函数 p 相对于参考测度 μ（通常是勒贝格测度）。考虑概率密度函数 p(·; θ) 依赖于一个有限维参数集 θ \in Θ \subset RK，其中 Θ 表示 θ 的所有可能值，K \in N。统计问题是在给定 y 的情况下恢复未知的 θ = (y1, . . . , yn)>，其中 yi 是 p(·; θ) 的独立副本。在观察到的数据样本 y 处评估 p(·; θ) 会给出一个实值函数 p(y; θ)，称为似然函数。在大多数情况下，这种似然函数是明确已知的，并且可以通过分析或数值方式对 (y; θ) 对进行评估。为了估计 θ，人们通常使用最大似然原理，该原理假设最合理的值是观察样本密度最大的值。使对数似然函数最大化的值 ^ θ= ^ θ(y) 称为最大似然估计 (MLE)：

θ=argmaxθ∈Θp(y;θ)。^ θ = arg max θ\in Θ p(y; θ)。

然而，通常没有已知或可用的最大化问题的封闭形式解决方案，并且 MLE 只能通过数值优化找到。在更极端的情况下，例如依赖关系的许多非高斯模型，写下完整的可能性是不可能的，即使在中等维度上，高级数值优化策略也不切实际。我们建议使用深度神经网络来解决这个问题，深度神经网络是一种机器学习，它利用连接的多层模型集来学习复杂的输入-输出映射。这个想法是使用参数样本和数据对来学习函数 p(·; θ)。换句话说，我们的目标是找到 p(·; θ) 的有用近似值，它是 θ 和 y 之间预测关系的基础。我们的估计技术是一个深度神经网络，将 y 作为独立输入数据，将 θ 作为相关输出数据

F:y7→θ;θ=argmind{y,Fw(θ)}F : y 7→ θ; ^ θ = arg min d\{y, Fw(θ)\}

其中 Fw 是一个函数，它依赖于一组 l1 个深度 NN 权重和 l2 个 NN 偏差，换句话说，w = (w1, . . . wl1, b1, . . ., bl2)>。这些权重由训练决定，而 (2) 中的 d 表示 y 和 Fw(θ) 之间的距离度量。函数 Fw 类似于 arg max 操作，但是它通过深度神经网络来近似 y 和 θ 之间的关系，而不是最大化似然函数。在我们的方法中，我们寻找函数 F 来预测 θ 给定输入 y 的值。这个问题需要一个选择F的标准，即惩罚预测错误的损失函数。对于回归问题，均方误差 (MSE) 损失是最常见和方便的选择 (Friedman et al., 2001)，由下式给出

MSE(w)=E{θ−Fw(θ)}2MSE(w) = E\{θ − Fw(θ)\}2

优化算法使用梯度，通常通过自动微分计算（Paszke 等人，2017），将（3）中的损失函数应用于深度神经网络和数据的输出。这个想法是，深度 NN 找到一组权重 w，使得拟合根据输出的实现尽可能接近观察点，其中接近度由 MSE(w) 测量。对于线性模型，存在 (3) 的最小化问题的简单封闭形式解；否则，解决方案通常需要迭代方法。实施了围绕批量梯度下降法构建的几个数值优化器来执行此任务（例如，参见 Kingma 和 Ba (2014)）。

优化 Fw 的问题涉及为深度神经网络生成训练数据。为此，我们生成参数和数据对的训练样本 (θtrain j , ytrain j )J j=1，其中 θtrain j = (θ1,j, . . . θK,j)> 和 Ytrain j = (y1,j, . . yn,j)>，其中 θtrain j 对应于覆盖感兴趣的参数域 Θ 的配置是至关重要的。可以使用不同的方法来找到这样的配置。一种可能的方法是从观察到的数据 y 中获得 θ 的粗略估计（例如，使用近似最大似然法）并在这些数据的足够大的邻域中模拟 (θtrain j )J j=1。然后，在定义了参数空间 Θ 并生成了 (θtrain j )J j=1 的样本之后，这些可以用来模拟相应的数据样本 (ytrain j )J j=1。训练样本的数量取决于问题的难度，并直接影响估计器的质量（Juba 和 Le，2019）。然而，我们发现输入和输出的适当转换大大减少了获得良好精度所需的样本数量。一旦深度神经网络被训练好，它就可以用来完成我们执行统计推理的目标，即根据新的观察数据返回概率模型的参数。

2.2 卷积神经网络简介

在这里，我们详细介绍了 (2) 中 Fw 假设的具体形式，负责学习模拟数据（输入）和深度神经网络的统计模型参数（输出）之间的映射。深度神经网络的选择可能取决于输入和输出的类型。在以下部分中，我们分析空间二维数据；因此，我们专注于 CNN，这是一类最常应用于视觉图像的深度神经网络，它解释了信号中的局部依赖性。我们认为 Fw 是 m 个非线性函数的组合，使得 Fw(θ) = (fm ◦ fm−1 ◦ . . . ◦ f1)(θ)，其中每个 f (·) 是至少一个隐藏层的函数. CNN 的一般结构交替使用卷积层和池化层，最后一层完全连接（密集）。卷积运算能够对输入执行加权平均，以便网络学习在检测到输入中的特定空间模式时激活的过滤器。形式上，给定输入图像 y 和内核 k，离散卷积运算符由下式给出

K\[s1, s2] ∗ y\[s1, s2] = ∞ ∑ i=−∞ ∞ ∑ j=−∞ K\[i, j]y\[s1 + i, s2 + j],

其中 s1 和 s2 是图像的像素。在实践中，总和是有限的，因为内核 K 具有紧凑的支持并且取决于像素的数量。每个内核 K 是一个可训练权重矩阵（也称为滤波器），应用于输入图像的小区域，并且相同的内核在图像的各个维度上进行转换。当在卷积层中使用多个过滤器时，可以应用几个独立训练的内核。根据内核权重的值，可以获得与图像相关的不同属性，例如对象的边缘。我们指的是 Pinaya 等人。 (2020) 有关 CNN 的更多信息，并访问 https://tensorflow.org/api_docs 以获取有关卷积层的 TensorFlow 实现的说明。

3 最大稳定过程的概念证明

最大稳定分布通常用于研究时空记录的极端事件。作为概念证明，我们通过估计最大稳定过程的参数来举例说明第 2 节中的方法。这些过程以具有即使在中等维度上也难以计算的完全似然而闻名。但是，模拟是可管理的。通过这个例子，我们的目的是证明我们的方法可以通过从原始数据中学习映射来准确地恢复具有难以处理的可能性的模型的参数。我们在这里使用 Schlather (2002) 提供的算法来模拟在 SpatialExtremes R 包中实现的最大稳定模型。

在第 3.1 节中，我们简要概述了最大稳定分布和过程，以及在此类模型中估计参数的最常用方法。第 3.2 节和第 3.3 节分别概述了参数估计设置和估计过程中使用的 CNN。在第 3.4 和 3.5 节中，通过使用两个流行的最大稳定模型的模拟研究，评估了 CNN 估计器的性能并将其与基准进行比较。

3.1 定义和成对似然法

考虑一系列独立同分布 (i.i.d.) 的随机过程 Y1(s), Y2(s), . . ., 其中 s \in S \subset Rd 是一个空间站点。最大稳定过程 Z 是归一化逐点块最大过程（块大小为 m）的极限过程，前提是函数序列 am(s) > 0 和 bm(s) 存在：

Z(s) = lim m→∞ am(s)−1 \[ max m \{Y1(s), . . . , Ym(s)\} − bm(s) ] , s \in Rd,

换句话说，Z(s) 的独立副本的逐点最大值保留在同一个族中，直到一个位置和尺度函数。 Z(s) 的每个边缘分布都是广义极值 (GEV) 分布，其位置、尺度和形状参数可能取决于空间位置（De Haan 和 Ferreira，2007 年）。

我们现在介绍一种表示，它将作为以下部分中构建参数化最大稳定过程的基础。令 {ξi}i≥1 为 (0, ∞) 上强度为 dΛ(ξ) = ξ−2dξ 的泊松过程的点，并令 W1(s), W2(s), . . .是非负随机过程的独立副本，使得 W (s) ≥ 0，均值等于 1。假设过程 Wi 和泊松过程的点 {ξi}i≥1 是独立的。然后

Z(s)=maxi≥1ξiWi(s),s∈S,Z(s) = max i≥1 ξiWi(s), s \in S,

是一个具有单位 Fr 的最大稳定过程：P{Z(s ≤ z)} = exp(−1/z), z > 0。W (·) 的合适选择会产生不同的最大稳定过程（Schlather , 2002).在接下来的部分中，我们将考虑两种不同的选择：Brown-Resnick (Kabluchko et al., 2009) 和 Schalther 的模型 (Schlather, 2002)。

• 当 (5) 中的 Wi(s) = exp{ i(s) − γ(s)} 时出现 Brown-Resnick 模型，其中 i(s) 是标准平稳高斯过程的独立副本，使得 (0) = 0 几乎肯定并且半变异函数 γ(h) = (|h|/λ)ν，其中 h 是空间间隔，λ > 0 是范围，而 ν \in (0, 2] 是平滑参数。Brown-Resnick过程在实践中很受欢迎，因为它们与其他选择相比具有灵活性，并且能够推广几个平稳的最大稳定模型，例如几何高斯模型（Davison 等人，2012）。有关其实际性能的精确证明，请参阅蒂博和奥皮茨 (2015)。

• 第二个模型是Schlather (2002) 提出的表征，其中Wi(s) = √2πmax{0, i(s)}，其中1(s), 2(s), . . .是具有单位方差和相关函数 ρ 的平稳高斯过程的独立副本，它是幂指数函数：ρ(h) = exp{−(|h|/λ)ν}, λ > 0, ν \in (0, 2] .

Schlather (2002) 表明，如果随机泊松过程 ξ 的支持（参见 (5)）包含在球 b(o; r) 中或者是由常数 C 统一界定。对于 Brown-Resnick 和 Schlather 的模型，这些条件不满足，因为高斯过程不是统一界定的。然而，Schlather (2002) 引入了常数 C 的近似值，这样如果 P[max{0; ξi(s)} > C] 足够小，模拟过程仍然准确。

从 (5) 可以看出，Z(s) 在有限的站点集合 {s1, . . . , sD} \subset S 由下式给出

p(Z(s1)≤z1,..,Z(sD)≤zD)=exp{−V(z1,...,zD)}p(Z(s1) ≤ z1, . . , Z(sD) ≤ zD) = exp\{−V (z1, . . . , zD)\}

其中指数函数 V (z1, . . . , zD) = E[max{W (s1)/z1, . . . , W (sD)/zD}], 满足同质性和边际约束 (De Haan et al., 1984)。通过微分分布 (6) 关于其变量 z1, . . . , zD, 可以推导出相应的密度函数，或完全似然

f(z1,...,zD)=exp{−V(z1,...,zD)}∑π∈PL∞l=1{−Vπl(z1,...,zD)},f (z1, . . . , zD) = exp\{−V (z1, . . . , zD)\} ∑ π\in P L ∞ l=1 \{−Vπl(z1, . . . , zD)\},

其中 P 表示所有分区的集合 π = {π1, . . . , πL} 的 {1, . . . , D} 和 Vπl 表示函数 V 相对于集合 πl 索引的变量的偏导数。等式 (7) 中的总和取自所有可能分区的集合，它等于 D 阶的贝尔数。即使对于适度的 D，这也会导致项的爆炸式增长，并且完全似然很快变得难以处理。卡斯特鲁乔等人。 (2016) 得出结论，使用现代计算机技术，完全似然推断仍然限于 D ≤ 12。然而，在简单的双变量情况下，密度在计算上是易于处理的，这解释了采用成对似然代替全似然的常见做法可能性（参见，例如，Padoan 等人 (2010)；Davis 等人 (2013)；Shang 等人 (2015)）。模型 (6) 的相应对数加权成对似然由下式给出

l(φ) = ∑ (j1,j2)\in P αj1,j2 \[ log\{V1(zj1, zj2)V2(zj1, zj2) − V12(zj1, zj2)\} − V1(zj1, zj2) ] ,

其中 φ \in Φ \subset Rp 是未知参数的向量，αj1,j2 ≥ 0 表示分配给 {j1, j2} 对的似然权重。为了简化符号，我们抑制了 V1、V2 和 V12 对 φ 的依赖性。权重通常仅针对特定距离内的成对位置设置为一个：h = |s1 − s2|< δ，其中 δ > 0 是选择的适当截止距离，例如，通过交叉验证（Ribatet , 2008).在第 j 个站点记录的块最大值用 zj 表示，其中 P={(j1,j2):1≤j1<j2≤D}P = \{(j1, j2) : 1 ≤ j1 < j2 ≤ D\}。对于 zj 的多个独立复制，(8) 中的最大成对似然估计量 φ 是强一致且渐近高斯分布的。这种错误指定的似然估计的渐近方差是三明治形式的（Padoan 等人，2010）。然而，计算这种三明治方差非常耗时，因为它需要一次处理四个点。也可以使用块自举技术评估 φ 的可变性。

3.2 建议的参数估计设置

为简单起见，我们在下面描述了我们提出的针对均匀边际模拟数据的估计方法，这样我们只估计依赖参数。我们模拟 I 独立重复 (θi, yi)I i=1（也称为测试数据），其中 θi = (λtest i , νtest i )> 是最大稳定过程的范围和平滑度，yi = {yi( s1), . . . , yi(sD)}> 是相应的实现（Brown-Resnick 或 Schalther 的）。我们在这里考虑 I = 16 个参数场景，其中 λtest i 和 νtest i 的值位于（不一定等距）网格上。对于每一对 (λtest i , νtest i )，在 [0, 20]2 中的 D = 252（即 25 × 25 图像）处生成最大稳定样本。然后我们使用两种不同的方法获得 θi 的估计值：

  1. CNN 估计器 ^ θCNN i = (^ λCNN i ,^ ν CNN i )>: 我们模拟训练数据 ytrain j 每个大小为 D = 252 的最大稳定模型，参数从 λtrain j∼ Unif(atrain λ , btrain λ ) 和 νtrain j ∼ Unif(atrain ν , btrain ν ), 对于 j = 1, . . . , 2000. atrain θ 和 btrain θ 的选择，对于 θ = λ, ν，由测试参数 θi 告知，以确保训练集覆盖感兴趣的区域。 CNN 使用 log(ytrain) 作为输入进行训练，其中 ytrain 是尺寸为 2000 × 25 × 25 的张量，{ log(λtrain), log ( ν train 2−ν train )}> 作为输出，尺寸为 2000 × 2. 此处使用对数作为方差稳定变换，以帮助解决训练期间的数值问题，并选择对应于 ν 的分母，使该参数不大于 2。训练完成后，CNN 用于根据对数变换的测试图像 yi 返回预测 ^ θCNN i，其中 i = 1，…。 . . , 16.

  2. Pairwise-likelihood estimator ^ θPL i = (^ λPL i ,^ ν PL i )>：通过 R 函数 fitmaxstab 使用成对似然（见（8））将模型拟合到数据。此函数使用另一个 R 函数 optim 执行优化，我们将方法设置为 L-BFGS-B。经过小型交叉验证研究后，我们发现除了增加拟合时间外，包括相距超过三个单位的对还会使估计变差。因此，我们只组合最多相隔三个单位且权重相等的对的似然贡献 αi1,i2 = 1，其余权重设置为零。我们通过为优化器提供围绕实际值的多个随机起始对来初始化参数。然后，我们从随机对中具有最高成对可能性的五对作为起点运行完整优化。我们使用具有成对似然最大化器的对作为初始值来显示结果。

表 1：CNN 模型总结。它是一个顺序模型，输入形状为 [–, 25, 25, 128] 并将其映射到形状为 [–, 2] 的两个标量值。

3.3 CNN架构

我们使用顺序 CNN，将模拟图像作为输入并将其映射到两个标量值。表 1 总结了用于使用来自 Brown-Resnick 和 Schlather 模型的模拟数据训练 CNN 的设置。我们报告层数以及它们如何转换输入张量的维度。例如，表1的第一行显示模型输入是形状为[–, 25, 25, 128]的张量，其中’–'表示任意数量的样本，数字’25’与维数有关图像，'128’代表过滤器的数量。然后，第一个卷积层使用 128 个内核大小为 3 × 3 的滤波器（见第二行）将输入张量转换为另一个输出形状为 [–, 13, 13, 128] 的张量。

拟议的 CNN 由三对具有不同数量过滤器的卷积组成。我们在网络的末端保留了三个密集层，分别有 8、16 和 2 个单元。 CNN 的可训练权重总数为 168,558。对于所有层，我们使用整流线性单元 (ReLU) 激活函数，这已被证明对回归类型模型有益（Hastie 等人，2009 年）。我们的算法的实施是通过使用 R 软件中内置的 TensorFlow/Keras 来执行的。 CNN权重是随机初始化的；为了训练网络，我们采用了广泛使用的 Adam 优化器（Kingma 和 Ba，2014 年），学习率为 0.01。训练进行了 32 个时期，在每个时期中，CNN 权重使用来自完整训练数据集的 40 个样本的批量大小进行更新，此处包含 J = 2000 个样本。

3.4 Brown-Resnick 模型的结果

在本节中，我们展示了 CNN 和成对似然估计 Brown-Resnick 模型参数的比较。第 3.1 节中给出了该模型的描述，第 3.2 节中提供了为这两种方法设置的参数估计的详细信息。测试集包含16对模拟参数和数据，其中参数为所有可能的组合(λtest, νtest)>，其中λtest\in [0.50,0.75,1.00,1.50]，νtest\in [0.80,1.05,1.30,1.55] .对于我们称之为场景的每个组合，我们模拟 50 i.i.d.各个 Brown-Resnick 过程的实现，将用于评估估计量的不确定性。涵盖模拟参数的实现示例如图 1 (a) 所示。

CNN 的训练集的输出值来自 λtrain j ∼ Unif(0.1, 3) 和 ν train j ∼ Unif(0.5, 1.9), j = 1, . . . , 2000, 以便它以合理的余量覆盖真实输出的范围。一旦网络收敛，我们就用它来预测测试数据上的 Brown-Resnick 参数。来自 CNN（绿色）和成对似然（红色）的估计变换范围 (log(λ)) 与变换平滑度 (log{ν/(2 − ν)})（原始比例轴上的数字）的散点图) 显示在图 2 (a) 中。每个图中的×符号代表真实情况。每个散点图包含 50 个独立估计的独立数据重复。在拟合成对似然时，我们注意到在多个起始值中使用最大化器的重要性（我们使用了 20 个围绕真值的值，并从具有最高成对似然的五对中进行了全面优化；有关详细信息，请参阅第 3.2 节）。该图显示 CNN 始终运行良好，而成对似然性能因场景而异。正如估计长期依赖性所预期的那样，该参数的较大值的可变性增加，特别是来自成对似然估计器。特别是，成对似然往往会低估平滑度参数。当范围和平滑度都很大时，这一点更为明显（参见图 2 (a) 的右下角）。对于较大的平滑度，有一些数据重复，成对似然无法估计两个参数（参见图 2 (a) 的最后一行）。对于具有中等范围的较粗糙的字段，估计值存在很大的可变性，在大多数情况下会低估平滑度。基于这项研究，CNN 导致总体估计偏差和方差比成对似然更小，大多数估计更接近实际参数值。

表 2 报告了使用 CNN 和成对似然法预测 Brown-Resnick 参数 (λ; ν) 的均方根误差 (RMSE)、平均绝对误差 (MAE) 和平均偏差值。每个分数都是 50 次重复和 16 种场景组合的结果。 Brown-Resnick 模型（第 2 列和第 3 列）的结果证实 CNN 估计器总体上优于成对似然。这对于平滑度参数更为明显，对于成对似然，偏差更高，RMSE 和 MAE 大约大 1.6 倍和 2.7 倍。成对似然估计产生的偏差与图 2 (a) 中的散点图一致，可能是因为仅使用成对变量的效率低下。此外，我们为每个测试数据考虑一个重复，增加这个数字可以使成对的估计更稳定。然而，研究高阶复合似然法超出了本文的范围，我们建议读者参考 Genton 等人。 (2011) 有关此方法的详细信息。

3.5 Schlather 模型的结果

与上一节中描述的 Brown-Resnick 模型的实验类似，我们模拟了 Schlather 模型的测试集，该模型包含 [0, 20] 中大小为 D = 252 的 16 个参数场景的 50 个独立副本。这里，用于模拟数据的参数 (λtest, νtest)> 由来自 λtest \in [0.50, 1.50, 2.0, 2.5] 和 νtest \in [0.80, 1.05, 1.30, 1.55] 的所有对组合组成。用于训练 CNN 的数据集包含相同大小的模拟图像，其中根据 λtrain j ∼ Unif(0.1, 3) 和 νtrain j ∼ Unif(0.5, 1.8) 生成范围和平滑度，j = 1, . . . , 2000. 图 1 (b) 显示了 Schlather 模型的模拟示例。与上一节中介绍的 Brown-Resnick 过程相比（参见图 1 (a)），此模型生成的数据噪声更大，呈现出更具挑战性的场景。

图 2 (b) 中的每个散点图显示了 50 次重复的 log(λ) 与 log{ν/(2 − ν)} 的估计值以及使用 CNN（绿色）或成对似然进行估计的不同场景（红色的）。与 Brown-Resnick 过程的分析相比（参见图 2 (a)），这两种方法都产生了变异性较小的估计，并且成对的偏差在这里不那么突出。例外情况是范围较小的情况，其中成对似然低估了一些数据集中的范围（参见图 2 (b) 的第一列）。 CNN 在预测参数和所有场景时都会产生小的偏差，并且对于大范围的情况，变异性略大。结果表明，CNN 估计器总体上表现良好，并不断改进成对似然结果。我们注意到 CNN 的准确性与训练值的范围无关，因为 CNN 下的所有估计都不会接近模拟参数值范围的边界。虽然在估计大多数参数值的范围时，成对似然只会产生很小的偏差，但在一些情况下，该方法表现不佳，产生的估计值比真实值小两到三倍。相比之下，CNN 可以很好地估计更大的范围，并且极不可能出现异常值。

表 2 证实了 CNN 相对于成对似然估计器的改进模式，其中 RMSE 和 MAE 值小了大约 60% 到 80%。成对似然击败 CNN 的唯一情况是平滑参数的平均偏差，这主要出现在小范围内，即当空间相关性相对较弱时（见图 2 (b)）。在这种情况下，CNN 对平滑度的估计似乎有一个不强烈依赖于实际平滑度的分布。计算是在 2.4 GHz 的 8 核英特尔酷睿 i9 处理器上进行的。从 CNN 的 BrownResnick 和 Schlather 模型中拟合和预测参数的耗用时间（以秒为单位）分别为 58 和 45，不包括生成训练样本的时间。相反，对于成对似然，这些时间分别为 1,305 和 14,960。基于 CNN 的参数估计，不包括生成训练样本的时间，将计算时间减少了 300 分之一。总的来说，与我们在 Brown-Resnick 模型中观察到的类似，CNN 提供了比成对似然更好的结果，需要大量更少的计算资源。 CNN的成功证实了我们之前的结论。它甚至更强烈地支持需要替代方法来克服当前用于拟合棘手模型的次优方法的困难。此外，这种结果有望在更高维度的设置中得到改善，其中成对似然估计器的效率损失更为显着。

表 2：使用 CNN 和成对似然法的 Brown-Resnick 和 Schlather 模型的 RMSE、MAE 和平均偏差。前三行每列中的两个数字分别代表估计范围和平滑度的分数。

4 温度数据的应用

略，详情参见原文。

5 讨论

在这项工作中，我们提出了一种新方法来估计基于深度神经网络的统计模型中的参数。作为概念证明，我们已经基于极具挑战性的最大稳定分布和过程推理问题测试了该方法。我们通过模拟和真实世界的数据说明了使用深度 NN 模型优于经典成对似然法的好处。我们对流行的 Brown-Resnick 和 Schlather 模型的模拟研究结果表明，与当前的计算加速方法相比，我们的方法具有更好的性能（偏差和方差可忽略不计）。与成对似然不同，所提出的深度学习方法依赖于完整模型的模拟，不涉及空间依赖性的近似。

我们遵循的范例是，如果可以从模型进行模拟，就可以进行推理。我们的方法使用模拟数据作为深度神经网络的输入，并学习预测统计参数。与众所周知的 ABC 类似，其思想是深度神经网络成为压缩数据和识别模型参数的替身。然而，与 ABC 的一个关键区别在于我们的方法自动将数据和参数映射到深度神经网络内部，从而消除了选择汇总统计数据和指标来比较数据的问题。另一个区别是，在 ABC 中，整个估计过程必须从头开始为每个单独的数据集重新运行，而在提议的框架中，估计是在单个训练阶段预先执行的，然后是能够处理多个测试集的更便宜的预测阶段。此外，我们的新方法受益于传统深度神经网络算法的最新进展，它可以从高度复杂的图像中自动学习特征，同时快速且易于应用。当模型中有许多参数时，我们的方法就会出现问题，直到为高维参数空间生成训练数据变得实际上难以处理。

在这项工作中，我们专注于对多变量极端情况进行建模，并遵循概念验证方法。我们可以使用此处介绍的相同想法来推断其他统计模型中的参数。特别是，依赖性的非高斯模型往往难以使用经典方法进行估计。未来工作的扩展包括非齐次泊松过程（Moller 和 Waagepetersen，2003 年）、流行病学模型（Lawson，2008 年）、离散随机人口动力学模型（Wood，2010 年）和随机集模型。

总而言之，我们表明深度神经网络模型易于构建、拟合速度快，并且在执行统计推理方面显示出有希望的结果，尤其是在常用技术失败时。

  • [1] Castruccio, S., Huser, R., and Genton, M. G. (2016). High-order composite likelihood inference for max-stable distributions and processes. Journal of Computational and Graphical Statistics, 25(4):1212–1229.
  • [2] Cooley, D., Naveau, P., and Poncet, P. (2006). Variograms for spatial max-stable random fields. In Dependence in Probability and Statistics, pages 373–390. Springer.
  • [3] Creel, M. (2017). Neural nets for indirect inference. Econometrics and Statistics, 2:36–49.
  • [4] Cremanns, K. and Roos, D. (2017). Deep Gaussian covariance network. arXiv preprint arXiv:1710.06202.
  • [5] Davis, R. A., Kl ̈ uppelberg, C., and Steinkohl, C. (2013). Statistical inference for max-stable processes in space and time. Journal of the Royal Statistical Society: SERIES B: Statistical Methodology, pages 791–819.
  • [6] Davison, A. C., Padoan, S. A., Ribatet, M., et al. (2012). Statistical modeling of spatial extremes. Statistical Science, 27(2):161–186.
  • [7] De Haan, L. et al. (1984). A spectral representation for max-stable processes. The Annals of Probability, 12(4):1194–1204.
  • [8] De Haan, L. and Ferreira, A. (2007). Extreme Value Theory: An Introduction. Springer Science & Business Media.
  • [9] Erhardt, R. J. and Smith, R. L. (2012). Approximate Bayesian computing for spatial extremes. Computational Statistics & Data Analysis, 56(6):1468–1481.
  • [10] Friedman, J., Hastie, T., Tibshirani, R., et al. (2001). The Elements of Statistical Learning, volume 1. Springer Series in Statistics.
  • [11] Genton, M. G., Ma, Y., and Sang, H. (2011). On the likelihood function of Gaussian max-stable processes. Biometrika, pages 481–488.
  • [12] Gerber, F. and Nychka, D. W. (2020). Fast covariance parameter estimation of spatial Gaussian process models using neural networks. Stat, page e382.
  • [13] Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Science & Business Media.
  • [14] Huang, W. K., Stein, M. L., McInerney, D. J., Sun, S., and Moyer, E. J. (2016). Estimating changes in temperature extremes from millennial-scale climate simulations using generalized extreme value (gev) distributions. Advances in Statistical Climatology, Meteorology and Oceanography, 2(1):79–103.
  • [15] Huser, R., Davison, A. C., and Genton, M. G. (2016). Likelihood estimators for multivariate extremes. Extremes, 19(1):79–103.
  • [16] Huser, R., Dombry, C., Ribatet, M., and Genton, M. G. (2019). Full likelihood inference for max-stable data. Stat, 8(1):e218.
  • [17] Jiang, B., Wu, T.-y., Zheng, C., and Wong, W. H. (2017). Learning summary statistic for approximate Bayesian computation via deep neural network. Statistica Sinica, pages 15951618.
  • [18] '
  • [19] Juba, B. and Le, H. S. (2019). Precision-recall versus accuracy and the role of large data sets. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4039–4048.
  • [20] Kabluchko, Z., Schlather, M., De Haan, L., et al. (2009). Stationary max-stable fields associated to negative definite functions. The Annals of Probability, 37(5):2042–2065.
  • [21] Kharin, V. V. and Zwiers, F. W. (2000). Changes in the extremes in an ensemble of transient climate simulations with a coupled atmosphere–ocean gcm. Journal of Climate, 13(21):37603788.
  • [22] Kingma, D. P. and Ba, J. (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  • [23] Lawson, A. B. (2008). Bayesian disease mapping: hierarchical modeling in spatial epidemiology. Chapman and Hall/CRC.
  • [24] Liu, K., Ok, K., Vega-Brown, W., and Roy, N. (2018). Deep inference for covariance estimation: Learning Gaussian noise models for state estimation. In 2018 IEEE International Conference on Robotics and Automation (ICRA), pages 1436–1443. IEEE.
  • [25] Mitchell, K. E., Lohmann, D., Houser, P. R., Wood, E. F., Schaake, J. C., Robock, A., Cosgrove, B. A., Sheffield, J., Duan, Q., Luo, L., et al. (2004). The multi-institution north american land data assimilation system (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modeling system. Journal of Geophysical Research: Atmospheres, 109(D7).
  • [26] Moller, J. and Waagepetersen, R. P. (2003). Statistical inference and simulation for spatial point processes. CRC Press.
  • [27] Padoan, S. A., Ribatet, M., and Sisson, S. A. (2010). Likelihood-based inference for max-stable processes. Journal of the American Statistical Association, 105(489):263–277.
  • [28] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., and Lerer, A. (2017). Automatic differentiation in pytorch.
  • [29] Pinaya, W. H. L., Vieira, S., Garcia-Dias, R., and Mechelli, A. (2020). Convolutional neural networks. In Machine Learning, pages 173–191. Elsevier.
  • [30] Radev, S. T., Mertens, U. K., Voss, A., Ardizzone, L., and K ̈ othe, U. (2020). Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE Transactions on Neural Networks and Learning Systems.
  • [31] Ribatet, M. (2008). SpatialExtremes: An R-package for modelling spatial extremes. R Package Version, pages 2–0.
  • [32] Ribatet, M., Cooley, D., and Davison, A. C. (2012). Bayesian inference from composite likelihoods, with an application to spatial extremes. Statistica Sinica, pages 813–845.
  • [33] Rudi, J., Bessac, J., and Lenzi, A. (2021). Parameter estimation with dense and convolutional neural networks applied to the Fitzhugh-Nagumo ODE. Mathematical and Scientific Machine Learning.
  • [34] Schlather, M. (2002). Models for stationary max-stable random fields. Extremes, 5(1):33–44.
  • [35] Shang, H., Yan, J., and Zhang, X. (2015). A two-step approach to model precipitation extremes in California based on max-stable and marginal point processes. The Annals of Applied Statistics, pages 452–473.
  • [36] Stephenson, A. and Tawn, J. (2005). Exploiting occurrence times in likelihood inference for componentwise maxima. Biometrika, 92(1):213–227.
  • [37] Tawn, J. A. (1988). Bivariate extreme value theory: models and estimation. Biometrika, 75(3):397–415.
  • [38] Thibaud, E., Aalto, J., Cooley, D. S., Davison, A. C., Heikkinen, J., et al. (2016). Bayesian inference for the Brown–Resnick process, with an application to extreme low temperatures. The Annals of Applied Statistics, 10(4):2303–2324.
  • [39] Thibaud, E. and Opitz, T. (2015). Efficient inference and simulation for elliptical Pareto processes. Biometrika, 102(4):855–870.
  • [40] Varin, C., Reid, N., and Firth, D. (2011). An overview of composite likelihood methods. Statistica Sinica, pages 5–42.
  • [41] Wood, S. N. (2010). Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310):1102–1104.
  • [42] Xia, Y., Mitchell, K., Ek, M., Cosgrove, B., Sheffield, J., Luo, L., Alonge, C., Wei, H., Meng, J., Livneh, B., et al. (2012a). Continental-scale water and energy flux analysis and validation for north american land data assimilation system project phase 2 (NLDAS-2): 2. validation of model-simulated streamflow. Journal of Geophysical Research: Atmospheres, 117(D3).
  • [43] Xia, Y., Mitchell, K., Ek, M., Sheffield, J., Cosgrove, B., Wood, E., Luo, L., Alonge, C., Wei, H., Meng, J., et al. (2012b). Continental-scale water and energy flux analysis and validation for the north american land data assimilation system project phase 2 (NLDAS-2): 1. intercomparison and application of model products. Journal of Geophysical Research: Atmospheres, 117(D3).
  • [44] Xu, X. and Reid, N. (2011). On the robustness of maximum composite likelihood estimate. Journal of Statistical Planning and Inference, 141(9):3047–3054.