大型数据集的空间统计竞赛

〖摘要〗随着空间数据集变得越来越大和笨重，对空间模型的精确推断在计算上变得令人望而却步。已经提出了各种近似方法来减少计算负担。尽管存在对这些近似方法的综合评论，但对于一些选定的方法，它们的性能比较仅限于中小型数据集。为了实现包含尽可能多的方法的全面比较，我们组织了大型数据集空间统计竞赛。本次竞赛具有以下创新特点：1）我们使用 ExaGeoStat 软件生成合成数据集，生成的实现数在 10 万到 100 万之间； 2）我们系统地设计了数据生成模型来表示具有广泛统计特性的空间过程，适用于高斯和非高斯情况； 3) 竞赛任务包括估计和预测，并通过多个标准评估结果； 4）我们公开了所有数据集和竞赛结果，以作为其他近似方法的基准。在本文中，我们公开了所有比赛细节和结果以及对比赛结果的一些分析。

〖原文〗 Huang, H. et al. (2021) ‘Competition on Spatial Statistics for Large Datasets’, Journal of Agricultural, Biological and Environmental Statistics, 26(4), pp. 580–595. Available at: https://doi.org/10.1007/s13253-021-00457-z.

1 引言

随着更好的观测技术和先进计算设备的发展，获取大型空间数据集变得更加容易和普遍。因此，空间统计中的统计推断在计算上变得具有挑战性。几十年来，当精确计算不可行时，人们提出了各种近似方法来对大规模空间数据进行建模和分析。然而，在文献中，使用那些提出的近似方法的统计推断的性能通常只用中小型数据集进行评估，可以获得精确的解决方案。然而，对于大型现实世界数据集，精确计算不再可行。近似方法的推论通常根据经验或通过拟合模型的预测准确性进行验证。

受比较不同近似方法的统计和计算效率挑战的推动，引发了一些开创性的工作:

Englund (1990) ^[9] 进行了一项非常早期的研究，调查不同空间模型的推断性能。该研究在具有数万个数据点的两个区域中使用了 Walker Lake 数据集（Srivastava，1987），他们观察到来自不同空间模型的插值结果存在相当大的差异。
Bradley 等 (2016)^[6] 回顾了各种空间预测因子，包括确定性和随机方法，并将它们应用于 $CO_2$ 的卫星测量。研究了三种不同大小的数据集（最大的一个包含数万个观察值），并且对不同方法的评估依赖于预测误差。
Heaton 等 (2019)^[12] 主要关注高斯过程 (GP)，并提出了一项竞赛，在竞赛中，研究小组使用他们选择的高斯过程近似方法或其他无模型算法方法对模拟和真实世界的数据集进行预测。两个数据集都包含 $150,000$ 个实现。来自高斯过程的模拟数据的协方差是已知的，并向竞赛参与者公开；真实世界的数据包括由卫星测量的地表温度，其真正的潜在协方差是未知的。基于预测误差或预测分布检查性能。
Wikle 等 (2017) ^[26] 讨论了一个通用任务框架的设计来比较不同的方法。此外，他们开发了一个网站，让研究人员可以将他们的预测脚本上传到网站服务器以获取 NASA OCO-2 数据，然后相关的预测性能将发布在网站排行榜上。

受这些作品的启发，我们组织了一场竞赛，即大型数据集空间统计竞赛，以涉及最近的方法并克服以往研究中存在的弱点。与以往的比赛或对比作品相比，我们的比赛有以下主要特点：

通过使用 ExaGeoStat软件（Abdulah 等人，2018a）^[2]模拟大数据集，我们可以更好地理解不同方法的统计效率。
数据集是根据各种空间模型模拟的，包括高斯和非高斯过程模型。对于高斯过程模型，数据集是用 Matern 协方差函数模拟的一组选定参数，代表空间随机场的广泛统计特性。非高斯空间数据集由 Tukey g-and-h 随机场（Xu 和 Genton，2017）^[27]生成，它概括了高斯过程以解释偏度和重尾。
竞赛任务包括估计和预测。检查高斯过程近似方法推断的模型与事实的偏离程度是有见地的。一项分项竞赛侧重于评估模型规格错误，我们使用平均效率损失 (MLOE) 和均方误差的平均规格错误 (MMOM) 作为标准，这两个标准均由 Hong 等 (2021) ^[13]提出。。对于侧重于预测的分项竞赛，使用均方根误差（RMSE）来评估预测精度。
我们在公共存储库中共享了所有准备好的数据集和竞赛结果。未来针对大型空间数据集的近似方法可以使用这些数据集作为基准数据，并将其性能与现有方法进行比较。此外，我们还发布了 ExaGeoStat 精确计算的模型参数估计结果和预测结果，可作为精确推断结果参考。

大赛于 2020 年 11 月 23 日启动，吸引了全球 29 支研究团队报名。这些注册团队包括空间统计界的活跃研究人员。比赛于 2021 年 2 月 1 日结束，21 支队伍成功提交成绩；有关提交结果的团队的详细信息，请参见补充材料中的表 S1。

本文的其余部分组织如下：在第 2 节中，我们提供了比赛的概述。在第 3 节中，我们简要讨论了参赛队伍在比赛中使用的方法。在第 4 节中，我们通过一些分析展示了比赛结果的细节。在第 5 节中，我们总结并给出最后的评论。

2 竞赛信息

2.1 比赛准备中的 ExaGeoStat 软件

略

2.2 比赛使用的数据集

略

2.3 分项比赛详情

分项竞赛 1a：检查不同方法在中等大型数据集上正确推断高斯过程模型的能力。我们在每个数据集 G1 – G16 中选择了 $90,000$ 个实现，并要求参与团队估计四个参数 $σ^2$ 、 $β > 0$ 、 $ν > 0$ 和 $τ^2$ 。用于评估性能的指标是 MLOE 和 MMOM（Hong 等，2021) 跨越不同的数据集。 MLOE 表征了当使用近似模型而不是真实模型进行预测时，预测效率的平均损失。 MMOM 表征了在近似模型下计算时均方误差的平均误指定。 MLOE 和 MMOM 的详细信息在补充材料的 S3 节中给出。

分项竞赛 1b：评估在从高斯模型生成的中等大型数据集支撑下的空间预测性能。对于每个数据集 G1 – G16，我们为参赛团队提供了 $10,000$ 个新位置，并要求他们根据分项竞赛 1a 中提供的 $90,000$ 个实现来预测这些位置。 RMSE 用于评估预测准确性。

分项竞赛 2a：侧重于非高斯数据的预测，要求参赛团队预测超过 $10,000$ 个新位置，条件是每个数据集 NG1 和 NG2 有 $90,000$ 个实现。

分项竞赛 2b：对更大的数据集进行建模。选择了一个高斯数据集（数据集 G5）和一个非高斯数据集（数据集 NG1）。我们将条件数据规模增加到 $900,000$ ，参赛团队需要预测超过 $100,000$ 个新位置。

分项竞赛情况

2.4 评估

略

2.5 结果

每个提交的完整竞赛结果在补充材料的表 S2 中给出。使用了许多近似方法，我们在第 3 节中提供了它们的简要总结。为了突出表现最好的人，每个分项竞赛中的前三名提交（ 分项竞赛 2a 中的四个提交由于三个团队之间的平局或第二名）如下，得分最高的队伍排在第一位：

分项竞赛 1a : [1] SpatStat-Fans, [2] GpGp, [3] RESSTE(CL/krig)
分项竞赛 1b : [1] RESSTE(CL/krig), [2] HCHISS, [3] Chile-Team
分项竞赛 2a : [1] RESSTE(Tukey-g-h-trans-GPGP), [3] GpGp(quick), [3] HMatrix, [3] RESSTE (nonpara-trans-GPGP)
分项竞赛 2b : [2] RESSTE(nonpara-trans-GPGP), [2] RESSTE(Tukey-g-htrans-GPGP), [2] Tohoku-University

除了竞赛提交，我们还使用 ExaGeoStat 查看 分项竞赛 1a 和 1b 中精确计算的排名。 分项竞赛 1a 和 1b 的扩充顶级列表如表 4 所示，其中 分项竞赛 1a 中的 “ExaGeoStat(estimated-model)” 表示我们使用 ExaGeoStat 通过精确计算最大化完全似然来估计 Matern 协方差参数；在 分项竞赛 1b 中，这意味着我们使用关联的估计模型进行精确计算的预测。 分项竞赛 1b 中的 “ExaGeoStat(true-model)” 表示使用真实模型进行预测，ExaGeoStat 精确计算。

Table04

我们观察到 ExaGeoStat(estimated-model) 在 分项竞赛 1a 中设法找到了最接近真实情况的模型，但在 分项竞赛 1b 中预测性能略差于 RESSTE(CL/krig)。这表明使用更接近事实的模型并不能保证给定实现的点预测性能更准确。然而，ExaGeoStat(estimated-model) 应该在空间过程的多个实现中产生平均最好的预测。当使用真实参数值进行预测时（ ExaGeoStat(true-model) ），分数从 79 分下降到 72 分，是 分项竞赛 1b 中的最好成绩。

3 方法分类

在本节中，我们不打算对现有的近似方法进行详尽的文献综述，仅简要讨论参与者在比赛中使用的方法。详细的文献综述可以在 Sun 等 (2012) 和 Heaton 等（2019）中找到

3.1 主要方法

（1）组合似然法

此类方法将联合似然近似为一组分量似然的加权乘积（Varin 等人，2011 年）^[24]。例如，Vecchia 的近似框架使用一系列条件似然，其中条件集是稀疏选择的（Vecchia，1988）^[25]。成对似然法将每对观测值的似然性作为成分似然性 (Varin, 2008) ^[23]。因此，可以用较少的计算得到组合似然中的每个分量。

Teams Among-Stats、Chile-Team、ExtStat、GpGp、HCHISS、RESSTE 等团队提交了组合似然近似法的结果。

（2）低秩近似方法

低秩逼近方法一般将整个随机过程投影到某个低维空间，并使用低秩表示作为代理来逼近原始过程。例如，预测过程 (Banerjee 等, 2008) ^[5] 在空间域中放置节点，并且以节点实现为条件的原始过程的期望被用作替代。 Fixed rank kriging (Cressie and Johannesson, 2008) ^[7] 使用少量基函数来表示过程，从而可以通过对维数更小的矩阵求逆来获得精度矩阵。

ExtStat、UOW 等团队提交了低秩近似方法的结果。

（3）稀疏近似方法

协方差锥化 (Furrer 等, 2006 ^[10]; Kaufman 等, 2008 ^[16]) 将具有紧凑支持的相关函数与原始协方差函数相乘，使得远距离位置的相关性缩小为零，并在协方差矩阵中引入稀疏性。
对于精度矩阵，高斯马尔可夫随机场自然会在精度矩阵中产生稀疏结构 (Rue 等, 2009) ^[19]。

ExtStat 团队使用提交了这种技术的结果。

（4）其他方法

层次矩阵方法 (Litvinenko 等, 2019) ^[17] 将层次矩阵近似格式应用于协方差矩阵。然后，协方差矩阵的非对角块采用低秩表示，这使得协方差矩阵可以以较低的计算成本求逆。HMatrix 团队通过这种方式提交了结果。
协方差函数的全尺度近似（Sang 和 Huang，2012）^[20] 结合了协方差锥化和预测过程模型，同时考虑了小尺度和大尺度的空间依赖性。SpatStat-Fans 团队使用了这种方法。
多分辨率近似 (Katzfuss, 2017) ^[15] 是预测过程模型或全尺度近似模型的扩展，其中具有层次结构的基函数用于捕获不同尺度的空间依赖性。Colorado-School-of-Mines 和 GPvecchia 团队提交了使用这种方法的论文。
最近邻高斯过程 (Datta 等, 2016) ^[8]将 Vecchia 近似扩展到基于过程的模型，以便使用统一框架估计参数和进行预测。 ExtStat 和 NNGP 团队在他们的提交中应用了这种方法。

3.2 主要参赛队伍及其使用的方法

在这里，我们为比赛中的顶级团队提供简要说明和设置。

（1） 分项竞赛 1a 中的高斯过程模型推断问题

SpatStat-Fans 采用平滑的 全尺度近似方法，将整个域划分为 10×10 的规则矩形块，结点集合位于 $20×20$ 的网格上。
GpGp 对 $30, 000$ 个观测值进行二次采样，然后使用 R 包 GpGp 以 $30$ 个最近邻为条件的 条件分解近似法（Guinness 等人，2021 年）^[11]。
RESSTE(CL/krig) 使用 组合似然法 来找到最佳协方差参数估计。

（2）分项竞赛 1b 的高斯过程预测问题

RESSTE(CL/krig) 使用即插即用克里金预测器，通过组合似然来推断参数。
HCHISS 通过过 Vecchia 近似估计协方差参数，然后使用以 $1, 000$ 个最近邻为条件的克里金法。
Chile-Team 采用高斯条件成对似然估计协方差参数，然后使用以 $800$ 个最近邻为条件的克里金法，。

（3）分项竞赛 2a 的非高斯响应预测问题

RESSTE(Tukey-g-h-trans-GPGP) 和 RESSTE(nonpara-trans-GPGP) 在 分项竞赛 2a 和 分项竞赛 2b 中应用了 Tukey g and-h 变换和非参数变换，使得变换后的数据分别近似为高斯分布，然后使用 R 包 GpGp 进行高斯预测。
HMatrix 使用层次矩阵逼近协方差矩阵，精度为 $10^{−6}$ 。
GpGp(quick) 使用 R 包 GpGp 中的 “matern nonstat var” 协方差函数，其中使用 $50$ 个基函数表示空间变化的协方差函数，协方差参数由 $10, 000$ 个随机样本估计 $20$ 个条件邻居；然后，通过克里金法与 $30$ 个条件邻居进行预测。

（4）分项竞赛 2b 的大规模预测问题（高斯和非高斯）

- RESSTE(Tukey-g-h-trans-GPGP) 和 RESSTE(nonpara-trans-GPGP) 在 分项竞赛 2b 中应用了 Tukey g and-h 变换和非参数变换，使得变换后的数据分别近似为高斯分布，然后使用 R 包 GpGp 进行高斯预测。
Tohoku-University 使用了协方差锥化，其中高斯过程的 Matern 协方差函数应用了通过交叉验证估计的参数。

论文的讨论部分将提供每个分项竞赛中顶级团队使用方法的更多细节。

4 竞赛结果分析

在本节中，我们提供了有关比赛结果的更多详细信息。

图（2） 展示了分项竞赛 1a 中所有团队提交的参数估计以及 ExaGeoStat 的精确计算以供比较。我们重点介绍了 ExaGeoStat 的结果以及分项竞赛 1a（SpatStat-Fans、GpGp、RESSTE(CL/krig)）和 1b（RESSTE(CL/krig)、HCHISS、Chile-Team）中的前三名。请注意，提交的 RESSTE(CL/krig) 在两个分项竞赛中均名列前三。除了 HCHISS 之外的所有提交都成功地非常精确地估计了块金参数。我们观察到，当过程更平滑（更大的平滑参数）并且具有更强的依赖性（更大的有效变程）时，参数估计通常更困难。在这种情况下，来自不同提交的偏基台和变程参数估计差异最大。

为了进行比较，我们还展示了提交的 HCHISS 和 Chile-Team 的模型推断结果，它们分别在分项竞赛 1b 中排名第 2 和第 3。然而，我们注意到他们的模型估计不如 SpatStat-Fans、GpGp 和 RESSTE(CL/krig)。

尽管 ExaGeoStat 总体上具有最准确的估计，但我们注意到对于数据集 G15 和 G16，具有精确计算的 ExaGeoStat 往往会高估偏基台和变程。

SpatStat-Fans 和 RESSTE(CL/krig) 显示出与精确计算结果相似的模式，但 GpGp 获得的估计更准确，更接近事实。补充材料中的图 S1 说明了绝对 MLOE 和 MMOM，其中我们观察到 GpGp 对于数据集 G15 和 G16 确实具有更小的绝对 MLOE 和 MMOM。估计参数的似然值也可以用于比较。当插入来自提交的参数估计值时，我们使用 ExaGeoStat 来计算精确的对数似然。补充材料中的图 S2 描述了来自提交的对数似然减去具有真实参数的对数似然。对于那些具有较小对数似然的方法，例如 Chile-Team 和 HCHISS，这意味着它们由于近似而未能找到似然的最大值。那些具有较高值的，例如 ExaGeoStat 和 SpatStat-Fans，可能已经获得了给定数据集的最优估计。

图 2：子竞赛 1a 中所有团队的参数估计的箱线图（为清楚起见，未显示异常值）。突出显示了 ExaGeoStat、SpatStat-Fans、GpGp、RESSTE(CL/krig)、Chile-Team 和 HCHISS 的真实值和估计值。突出显示的提交（真相除外）的图例顺序遵循它们在子竞赛 1a 中的排名。数据集 G9 - G16 分别与 G1 - G8 共享相同的协方差结构，除了块块。

图（3） 显示了分项竞赛 1b 中每个数据集的不同提交的 RMSE。我们强调与之前讨论的分项竞赛 1a 相同的提交，包括分项竞赛 1b 和 1a 中的前三名提交。此外，我们使用 ExaGeoStat 在分项竞赛 1a 中使用真实参数和 ExaGeoStat 的估计进行精确计算进行预测，并给出并突出显示相应的 RMSE。在 图（3） 的顶部面板中，我们使用箱线图来总结不同数据集的整体预测性能。由于分项竞赛 1a 和 1b 中顶级团队的 RMSE 无法使用箱线图尺度很好地区分，我们还在 图（3） 的底部面板中用条形图显示了他们的 RMSE，以便更好地比较这些顶级团队。我们观察到，当块存在时，RMSE 通常较大，因为数据具有较高的噪声水平。值得注意的是，分项竞赛 1a 中表现最好的 SpatStat-Fans 和 GpGp 成功地找到了出色的参数估计。然而，与其他突出显示的提交相比，他们更好的推断模型并没有带来更好的整体预测。一个可能的原因是他们的近似在克里金法中不充分，即使他们使用的基础模型更准确。事实上，GpGp 只使用了 50 个最近的邻居作为每个预测的条件集，而 HCHISS 使用了 1000 个最近的点。这表明模型推断和考虑的邻居数量对于局部克里金预测都很重要；很难说邻居的数量在多大程度上很重要。

图 3：子竞赛 1b 中每个数据集中所有提交的 RMSE 箱线图。还给出了 ExaGeoStat 在子竞赛 1a 中使用真实参数和估计参数的 ExaGeoStat 预测。在图例中，突出显示的提交按照它们在子竞赛 1b 中的排名顺序列出。子竞赛 1a 和 1b 中顶级团队的 RMSE 突出显示并显示在条形图中。数据集 G9 - G16 分别与 G1 - G8 共享相同的协方差结构，除了块块。

分项竞赛 2a 和 2b 中的 RMSE 总结分别在补充材料的图 S3 和 S4 中给出，我们在其中突出显示了两个分项竞赛中的顶级团队。表现最好的包括将 Tukey g-and-h 变换和非参数变换应用于高斯过程，以及基于推断（非平稳）高斯过程模型的其他局部克里金预测。

5 讨论

在本次比赛中，我们创建并发布了一组不同设计的基准数据。我们知道用于生成数据集的真实参数以及 ExaGeoStat 的精确最大似然估计，可用于研究未来提出的方法。出于实际原因，我们在本次比赛中仅选择并使用了生成的高斯过程数据集的子集。具有一百万个空间位置的完整数据集已在 https://doi.org/10.25781/KAUST-8VP2V 上公开提供，以便于在未来的研究中使用。未来的近似方法可以使用此存储库作为工具，根据不同参赛团队的提交来评估它们的性能，并在本次比赛中使用 ExaGeoStat 进行精确推断（ExaGeoStat 在分项竞赛 1a 中的精确最大似然估计的详细摘要也是在补充材料的表 S3 中给出）

我们没有在本次比赛中比较计算时间，因为参赛团队在自己的机器上对数据进行建模，执行时间无法直接比较。但我们在补充材料的图 S5 中总结了所有提交的分项竞赛 1b、2a 和 2b 中进行预测的执行时间。以 $90,000$ 次观察为条件进行 $10,000$ 次预测的中位时间对于高斯数据（在分项竞赛 1b 中）约为 60 秒，对于非高斯数据（在分项竞赛 2a 中）约为 430 秒。对于分项竞赛 2b 中更大的数据集，以 $900,000$ 次观察为条件进行 $100,000$ 次预测的中位时间约为 $2,700$ 秒。

我们还注意到，理想情况下需要复制具有相同设置的数据集，以便从统计的角度更好地评估不同的方法。然而，本次比赛使用的数据集已经非常庞大，许多团队无法进行多次重复的推断和预测。为了让大多数参赛团队都能参加比赛，我们在每个设置中只使用了一个副本。尽管如此，我们考虑的各种协方差设置为大型空间数据建模提供了一个公平的比较。

几十年来，由于大型空间数据集的普遍存在带来的挑战，大空间数据问题一直是一个活跃的研究领域，这些数据集通常包含数百万个观测值，例如遥感气候数据或数值模型输出。 “大数据” 研究领域因实际应用中空间数据的规模而得到推进。除了为更大的空间数据集开发高效和准确的方法外，最近的研究还集中在多变量空间和时空数据上，这些数据的大小可以显著放大。然后，预测问题将包括单个或多个变量的空间插值和时间预测。提供一个统一的框架来理解现有近似方法的性能在模拟和评估中更具挑战性，但对于建议未来的研究方向至关重要。

参考文献

[1] Abdulah, S., Y. Li, J. Cao, H. Ltaief, D. E. Keyes, M. G. Genton, and Y. Sun (2019). ExaGeoStatR: A package for large-scale geostatistics in R. arXiv preprint arXiv:1908.06936 .
[2] Abdulah, S., H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes (2018a). ExaGeoStat: A high performance unified software for geostatistics on manycore systems. IEEE Transactions on Parallel and Distributed Systems 29 (12), 2771–2784.
[3] Abdulah, S., H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes (2018b). Parallel approximation of the maximum likelihood estimation for the prediction of large-scale geostatistics simulations. In 2018 IEEE International Conference on Cluster Computing (CLUSTER), pp. 98–108.
[4] Abdulah, S., H. Ltaief, Y. Sun, M. G. Genton, and D. E. Keyes (2019). Geostatistical modeling and prediction using mixed precision tile Cholesky factorization. In 2019 IEEE 26th International Conference on High Performance Computing, Data, and Analytics (HiPC), pp. 152–162.
[5] Banerjee, S., A. E. Gelfand, A. O. Finley, and H. Sang (2008). Gaussian predictive process models for large spatial data sets. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70 (4), 825–848.
[6] Bradley, J. R., N. Cressie, and T. Shi (2016). A comparison of spatial predictors when datasets could be very large. Statistics Surveys 10, 100–131. CHAMELEON (2021, January). The Chameleon project. Available at https://project. inria.fr/chameleon.
[7] Cressie, N. and G. Johannesson (2008). Fixed rank kriging for very large spatial data sets. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70 (1), 209–226.
[8] Datta, A., S. Banerjee, A. O. Finley, and A. E. Gelfand (2016). Hierarchical nearestneighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association 111 (514), 800–812.
[9] Englund, E. J. (1990). A variance of geostatisticians. Mathematical Geology 22 (4), 417–455.
[10] Furrer, R., M. G. Genton, and D. Nychka (2006). Covariance tapering for interpolation of large spatial datasets. Journal of Computational and Graphical Statistics 15 (3), 502–523
[11] Guinness, J., M. Katzfuss, and Y. Fahmy (2021). GpGp: Fast Gaussian Process Computation Using Vecchia’s Approximation. R package version 0.3.2.
[12] Heaton, M. J., A. Datta, A. O. Finley, R. Furrer, J. Guinness, R. Guhaniyogi, F. Gerber, R. B. Gramacy, D. Hammerling, M. Katzfuss, F. Lindgren, D. W. Nychka, F. Sun, and A. Zammit-Mangion (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics 24 (3), 398–425. HICMA (2021, January). The HiCMA project. Available at https://github.com/ecrc/ hicma.
[13] Hong, Y., S. Abdulah, M. G. Genton, and Y. Sun (2021). Efficiency assessment of approximated spatial predictions for large datasets. Spatial Statistics, to appear.
[14] Johnson, S. G. (2014). The NLopt nonlinear-optimization package. Available at https://github.com/stevengj/nlopt.
[15] Katzfuss, M. (2017). A multi-resolution approximation for massive spatial datasets. Journal of the American Statistical Association 112 (517), 201–214.
[16] Kaufman, C. G., M. J. Schervish, and D. W. Nychka (2008). Covariance tapering for likelihood-based estimation in large spatial data sets. Journal of the American Statistical Association 103 (484), 1545–1555.
[17] Litvinenko, A., Y. Sun, M. G. Genton, and D. E. Keyes (2019). Likelihood approximation with hierarchical matrices for large spatial datasets. Computational Statistics & Data Analysis 137, 115–132.
[18] R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
[19] Rue, H., S. Martino, and N. Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71 (2), 319–392.
[20] Sang, H. and J. Z. Huang (2012). A full scale approximation of covariance functions for large spatial data sets. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74 (1), 111–132.
[21] Srivastava, R. M. (1987). A non-ergodic framework for variograms and covariance functions. Master’s thesis, Stanford University, Stanford, CA.
[22] Sun, Y., B. Li, and M. G. Genton (2012). Geostatistics for large datasets. In E. Porcu, J.M. Montero, and M. Schlather (Eds.), Advances and Challenges in Space-time Modelling of Natural Events, Volume 207, Chapter 3, pp. 55–77. Springer.
[23] Varin, C. (2008). On composite marginal likelihoods. Advances in Statistical Analysis 92 (1), 1–28.
[24] Varin, C., N. Reid, and D. Firth (2011). An overview of composite likelihood methods. Statistica Sinica 21, 5–42.
[25] Vecchia, A. V. (1988). Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society: Series B (Methodological) 50 (2), 297–312.
[26] Wikle, C. K., N. Cressie, A. Zammit-Mangion, and C. Shumack (2017). A common task framework (ctf) for objective comparison of spatial prediction methodologies. Stats & Data Science Views. Available at https://www.statisticsviews.com/article/a-common-task-framework-ctf-for-objective-comparison-of-spatialprediction-methodologies.
[27] Xu, G. and M. G. Genton (2017). Tukey g-and-h random fields. Journal of the American Statistical Association 112 (519), 1236–1249.

附录

S1 补充表格

（1）参赛队伍名单

TableS1-1
TableS1-2
TableS1-3

（2）分项比赛排名

表 S2：每个子比赛中每个提交的排名。表中“S1a”、“S1b”、“S2a”、“S2b”分别表示分项比赛1a、1b、2a、2b。粗体值表示每个子比赛的获胜者。

TableS2

（3）ExaGeoStat 的精确最大似然估计

表 S3：ExaGeoStat 对子竞赛 1a 中每个数据集的精确最大似然估计及其相应的 MLOE、MMOM 和对数似然。

TableS3

S2 补充图件

Figure-s1

图 S1：子竞赛 1a 中所有团队的绝对 MLOE 和 MMOM。突出显示了 ExaGeoStat、SpatStat-Fans、GpGp、RESSTE(CL/krig)、Chile-Team 和 HCHISS 的结果。灰点是来自其他提交的指标。在图例中，突出显示的提交按照它们在子竞赛 1a 中的排名顺序列出。数据集 G9 G16 分别与 G1 – G8 共享相同的协方差结构，除了块块。

FigureS2

图 S2：ExaGeoStat、SpatStat-Fans、GpGp、RESSTE(CL/krig)、智利团队和 HCHISS 的对数似然减去子竞赛 1a 中真实参数的对数似然。在图例中，突出显示的提交按照它们在子竞赛 1a 中的排名顺序列出。数据集 G9 - G16 分别与 G1 - G8 共享相同的协方差结构，除了块块。

FigureS3