最近邻高斯过程模型

【摘 要】 用于分析地统计数据的空间过程模型需要进行计算，随着空间位置的数量变大，这些计算变得令人望而却步。本文开发了一类高度可扩展的最近邻高斯过程 (NNGP) 模型，为大型地统计数据集提供完全基于模型的推断。我们确定最近邻高斯过程是一个定义明确的空间过程，它提供合法的有限维高斯密度和稀疏精度矩阵。我们将最近邻高斯过程作为稀疏归纳先验嵌入到丰富的分层建模框架中，并概述了如何在不存储或分解大型矩阵的情况下执行计算高效的马尔可夫链蒙特卡罗 (MCMC) 算法。该算法每次迭代的浮点运算 (flops) 与空间位置的数量成线性关系，从而呈现出可观的可扩展性。我们使用模拟研究说明了最近邻高斯过程相对于竞争方法的计算和推断优势，并且还分析了美国森林资源清查数据集中的森林生物量，其规模超过了其他降维方法。本文的补充材料可在线获取。

【原 文】 Datta, A. et al. (2016) ‘Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets’, Journal of the American Statistical Association, 111(514), pp. 800–812. Available at: https://doi.org/10.1080/01621459.2015.1044091.

【阅后感】

1 引言

随着地理信息系统和用户友好软件功能的不断增强，当今统计学家经常会遇到地理参考数据集，其中包含大量不规则位置的多变量观测值。这反过来激发了人们对位置参考数据统计建模的极大兴趣；请参见 Stein (1999)^[39]、Moller 和 Waagepetersen (2003) ^[31]、Schabenberger 和 Gotway (2004) ^[35]、Cressie 和 Wikle (2011) ^[8] 以及 Banerjee、Carlin 和 Gelfand (2014) ^[1] 等书籍以了解各种方法和应用。

空间过程模型（Spatial process models ） 在兴趣区域 $\mathcal{D}$ 上使用隐随机场 $\{\mathbf{w}(\mathbf{s}) : s \in \mathcal{D}\}$ 引入观测值之间的空间依赖性，该随机场可以指定任意有限随机变量集合的联合概率分布。例如，零均值的高斯过程确保 $\mathbf{w} = (\mathbf{w}(\mathbf{s}_1), \mathbf{w}(\mathbf{s}_2), \ldots , \mathbf{w}(\mathbf{s}_n))^{\prime} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}(\boldsymbol{\theta}))$，其中 $\mathbf{C}(\boldsymbol{\theta})$ 代表由未知参数 $\boldsymbol{\theta}$ 索引的一个协方差矩阵族。这种空间过程提供了一个丰富的建模框架，并被研究人员广泛使用以帮助理解科学中的复杂空间现象。

然而，该模型框架的拟合通常涉及求解 $\mathbf{C}(\boldsymbol{\theta})$ 的逆和行列式，需要 $n^3$ 阶浮点运算（flops）和 $n^2$ 阶的存储。当 $n$ 很大并且 $\mathbf{C}(\boldsymbol{\theta})$ 没有可优化利用的特殊结构时，这变得令人望而却步。

从广义上讲，对大型空间数据集进行建模是从 “利用低秩模型” 或 “使用稀疏性” 开始的。

1.1 低秩方法

低秩模型 试图在低维子空间上构建空间过程，一般通过在较小的 $r \ll n$ 个位置（ “结点（knots）” 或 “中心”）集合上，对原始空间过程（简称为父过程）的实现进行回归。例如：

• Higdon 2001 ^[24]；
• Kammann 和 Wand 2003 ^[25]；
• Rasmussen 和 Williams 2005 ^[32]；
• Stein 2007 ^[41]、2008 ^[42]；
• Banerjee 等 2008 ^[3]；
• Crainiceanu 等 2008 ^[6]；
• Cressie 和Johannesson 2008 ^[7]；
• Finley、Banerjee 和 McRoberts 2009 ^[15]。

低秩模型拟合的算法成本会从 $\mathcal{O}(n^3)$ 减少到 $\mathcal{O}(nr^2 + r^3) \approx \mathcal{O}(nr^2)$ 。不过，当 $n$ 很大时，经验调查表明 $r$ 必须足够大才能充分地近似父过程，因此造成 $nr^2$ 也会变得过高（参见 第 5.1 节）。此外，当相邻观测值高度相关并且空间信号能够控制噪声时，低秩模型表现不佳 (Stein 2014) ^[43]。尽管经过 偏置调整（bias-adjusted） 的低秩模型通常会得到更好表现（Finley、Banerjee 和 McRoberts 2009 ^[15]；Banerjee 等 2010 ^[2]；Sang 和 Huang 2012 ^[34]），但它们往往会增加计算负担。

1.2 稀疏方法

稀疏方法的常见策略是 协方差矩阵锥化 和 精度矩阵的稀疏化。

• 协方差矩阵的锥化。它使用紧凑支撑的协方差函数，在协方差矩阵 $\mathbf{C}(\boldsymbol{\theta})$ 中引入稀疏性。此类方法对于模型的参数估计和响应变量的插值（“克里金法”）非常有效，但在对残差或隐过程的更一般性推断方面，尚未得到充分开发或探索。参见例如：

  • Furrer、Genton 和 Nychka，2006 年 ^[16]；
  • Kaufman、Scheverish 和 Nychka，2008 年 ^[27]；
  • Du、Zhang 和 Mandrekar，2009 年 ^[11]；
  • Shaby 和 Ruppert，2012 年 ^[36]。

• 精度矩阵的稀疏化。直接在精度矩阵（即协方差矩阵的逆矩阵， $\mathbf{C}^{-1}(\boldsymbol{\theta})$ ）中引入稀疏性，是另一种常用的稀疏化策略。此类方法大多涉及采用某种形式的简化似然来近似高斯过程的完全似然。与低秩过程不同，此类方法不一定能够扩展到任意位置的新随机变量。此类方法也可能没有对应的空间过程，这限制了对空间协方差参数的推断。此外，此类方法通过将估计输入从不同过程模型导出的插值器来实现任意位置的空间预测，这可能无法准确反映预测不确定性的估计。例如：

  • 马尔可夫随机场方法（例如 Rue 和 Held 2005^[33]; Lindgren 2011^[30] ）；
  • 条件分布分解方法（Vecchia 1988 ^[46]、1992 ^[47]；Stein、Chi 和 Welty 2004 ^[44]）；
  • 组合似然方法（例如，Bevilacqua 和 Gaetan 2014 ^[5]；Eidsvik 等 2014 ^[12]）。

1.3 本文贡献

我们的贡献是： 为空间过程（可能完全是无法观测的）的完全基于过程的推断，提供了可观的可扩展性。从 “有限维的稀疏似然”（指似然分解等方法） 过渡到 “稀疏归纳的空间过程”（指低秩方法） 可能非常复杂，我们的策略是：

• 首先，使用从有向无环图构造的特定邻居集，在有限维概率模型中引入了稀疏性；
• 然后，利用上述邻居集合，将原有限维概率模型扩展为一个在不可数集合上的有效空间过程。

我们称此过程为 最近邻高斯过程(Nearest-Neighbor Gaussian Process， NNGP)，其有限维实现具有封闭形式的稀疏精度矩阵。

在 Vecchia (1988) ^[46]、Stein、Chi 和 Welty (2004) ^[44]、Emory (2009) ^[13]、Gramacy 和 Apley (2014) ^[21]、Gramacy、Niemi 和 Weiss (2014) ^[23] 以及 Stroud、Stein 和Lysen (2014)^[45] 中，尽管稀疏性已经被有效地用于近似高成本似然，但迄今为止，尚没有一个完全基于过程的建模和推断框架，而 本文提出的最近邻高斯过程，将模型参数估计、响应变量预测、底层空间过程插值纳入了一个高度可扩展的统一框架，填补了这一空白，并丰富了现有方法的推断能力。

为了展示其完整推断能力，我们将最近邻高斯过程部署为贝叶斯框架中空间过程的稀疏归纳先验（即作为分层模型中的第二个阶段：过程模型）。与低秩过程方法不同，最近邻高斯过程始终指定非退化的有限维分布，并且适用于空间随机过程的任意类型分布。因此，它可以被用来模拟未被实际观测到的隐过程。这种建模方式可以在第二阶段 （即分层建模框架中的过程模型阶段）为随机效应（例如截距或系数）提供结构化的依赖性，并且不要求第一阶段（即分层建模框架中的数据模型阶段）必须是高斯的。

我们在空间变系数回归框架（Gelfand 等，2003 年 ^[20]；Banerjee 等，2008 年 ^[3]）中使用了一个多元最近邻高斯过程，可以很方便地获得所有模型参数以及空间过程的完整后验。

最后，我们在一个规模较大的林业案例中（其规模大到即便低秩过程在高性能计算环境中也无法实现），展示了最近邻高斯过程为空间变系数回归模型提供的过程推断实例。

1.4 文章结构

• 第 2 节使用多元高斯过程制定最近邻高斯过程。

• 第 3 节概述了非常灵活的分层建模设置中的贝叶斯估计和预测。

• 第 4 节讨论替代最近邻高斯过程模型和算法。

• 第 5 节介绍了模拟研究，以突出最近邻高斯过程的推断优势，还分析了来自美国农业部庞大数据集的森林生物量。

• 第 6 节简短地总结了本文和对未来工作的启示。

2 最近邻高斯过程

2.1 稀疏有向无环图上的高斯密度

我们将考虑 $\mathbb{R}^d$ 上的一个 “$q$ -变量” 空间过程（注：指具有 $q$ 个响应变量或对应的随机过程）。

设 $\mathbf{w}(\mathbf{s}) \sim \mathcal{GP}(\mathbf{0}, \mathbf{C}(·, · | \boldsymbol{\theta}))$ 为一个零均值的 “$q$ -变量” 高斯过程，其中对于所有的 $\mathbf{s} \in \mathcal{D} \subseteq \mathbb{R}^d$，有 $\mathbf{w}(\mathbf{s}) \in \mathbb{R}^q$ 。该过程完全由有效的 互协方差函数（cross-covariance function） $\mathbf{C}(·, · | \boldsymbol{\theta})$ 指定，它将 $\mathcal{D} × \mathcal{D}$ 中的一对位置 $\mathbf{s}$ 和 $\mathbf{t}$ 映射到一个 $q × q$ 实值矩阵 $\mathbf{C}(\mathbf{s,t})$ 中，其中元素值为 $\text{cov}\{w_i(\mathbf{s}), w_j(\mathbf{t})\}$。这里，$\boldsymbol{\theta}$ 表示与互协方差函数相关的参数。设 $\mathscr{S}= \{ \mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_k\}$ 是 $\mathcal{D}$ 中不同位置组成的一个固定集合，我们称之为 参考集（reference set）。所以，$\mathbf{w}_{\mathscr{S}} \sim \mathcal{N}(\mathbf{0}, \mathbf{C}_{\mathscr{S}}(\boldsymbol{\theta}))$，其中 $\mathbf{w}_{\mathscr{S}} = (\mathbf{w}(\mathbf{s}_1)^{\prime}, \mathbf{w}(\mathbf{s}_2)^{\prime}, \ldots , \mathbf{w}(\mathbf{s}_k)^{\prime})^{\prime}$，并且 $\mathbf{C}_{\mathscr{S}}(\boldsymbol{\theta})$ 是一个大小为 $qk× qk$ 的正定分块矩阵，其中共有 $k \times k$ 个块，每个块的大小为 $q \times q$ 对应于 $\mathbf{C}(\mathbf{s}_i, \mathbf{s}_j)$。为便于表示，本文后面将 $\mathbf{C}_{\mathscr{S}}(\boldsymbol{\theta})$ 简写为 $\mathbf{C}_{\mathscr{S}}$ ，对 $\boldsymbol{\theta}$ 的依赖被隐含了，以后的所有空间协方差矩阵都将采用类似符号表示。

参考集 $\mathscr{S}$ 不需要与观测位置集合重合或成为观测位置集的一部分，因此 $k$ 不一定等于 $n$，尽管我们稍后表明观测位置集可能是 $\mathscr{S}$ 方便而实际的选择。当 $k$ 很大时，参数估计在计算上变得很麻烦，甚至不可行，因为它需要求 $C_{\mathscr{S}}$ 的逆和行列式。我们的方法受益于将 $\mathbf{w}_{\mathscr{S}}$ 的联合密度表示为若干条件密度的乘积：

$$p(\mathbf{w}_{\mathscr{S}}) = p(\mathbf{w}(\mathbf{s}_1)) \cdot p(\mathbf{w}(\mathbf{s}_2) \mid \mathbf{w}(\mathbf{s}_1)) \ldots \cdot p(\mathbf{w}(\mathbf{s}_k) \mid \mathbf{w}(\mathbf{s}_{k − 1}), \ldots , \mathbf{w}(\mathbf{s}_1)) \tag{1}$$

并将 式 (1) 右侧较大的条件集合替换为经过精心挑选的、大小不超过 $m$ 的条件集合（ $m \ll k$ ）。 参见例如，Vecchia 1988 ^[46]；Stein、Chi 和 Welty 2004 ^[44]； Gramacy 和 Apley 2014 ^[21]；Gramacy、Niemi 和 Weiss 2014 ^[23]）。因此，对于每个 $\mathbf{s}_i \in \mathscr{S}$，我们都会对应一个较小的条件集 $N(\mathbf{s}_i) \subset \mathscr{S} -\{\mathbf{s}_i\}$，进而可以将上式重写为：

$$\tilde{p}(\mathbf{w}_{\mathscr{S}}) = \prod^{k}_{i=1} p(\mathbf{w}(\mathbf{s}_i) \mid \mathbf{w}_{N(\mathbf{s}_i)}) \tag{2}$$

其中 $\mathbf{w}_{N(\mathbf{s}_i)}$ 是将 $\mathbf{w}(\mathbf{s})$ 在 $N(\mathbf{s}_i)$ 上的实现堆叠而成的向量。

设 $N_{\mathscr{S}} = \{ N(\mathbf{s}_i); i= 1, 2, \ldots , k\}$ 是 $\mathscr{S}$ 上所有条件集组成的集合。我们可以将 $\{ \mathscr{S} , N_{\mathscr{S}} \}$ 视为一个有向图 $\mathscr{G}$，其中 $\mathscr{S}= \{\mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_k \}$ 是结点的集合，$N_{\mathscr{S}}$ 是有向边的集合。对于任意两个结点 $\mathbf{s}_i$ 和 $\mathbf{s}_j$，如果存在来自 $\mathbf{s}_i$ 到 $\mathbf{s}_j$ 的有向边，则我们称 $\mathbf{s}_j$ 是 $\mathbf{s}_i$ 的 有向邻居（directed neighbor）。因此，$N(\mathbf{s}_i)$ 表示了 $\mathbf{s}_i$ 的有向邻居集合，以后我们称其为 $\mathbf{s}_i$ 的 “邻居集（neighbor set）”。有向图中的 “有向循环（directed cycle）” 指由结点 $\mathbf{s}_{i1}, \mathbf{s}_{i2}, \ldots , \mathbf{s}_{ib}$ 依次有向连接构成的链，其中链首和链尾为同一结点，即 $\mathbf{s}_{i1} = \mathbf{s}_{ib}$。不存在有向环的有向图被称为 “有向无环图（directed acyclic graph, DAG）”。

如果 $\mathscr{G}$ 是一个有向无环图，则根据如上定义，$\tilde{p}(\mathbf{w}_\mathscr{S})$ 是一个多元联合密度（参见在线附录 A1 或 Lauritzen (1996) ^[28] 以获得类似结论）。从联合多元密度 $p(\mathbf{w}_{\mathscr{S}})$ 开始，我们使用有向无环图 $\mathscr{G}$ 导出了新密度 $\tilde{p}(\mathbf{w}_{\mathscr{S}})$。虽然这适用于任何原始密度 $p(\mathbf{w}_{\mathscr{S}})$，但它在我们的上下文中特别有用，因为 $p(\mathbf{w}_{\mathscr{S}})$ 是多元高斯密度并且 $\mathscr{G}$ 足够稀疏。准确地说，令 $\mathbf{C}_{N(\mathbf{s}_i)}$ 为 $\mathbf{w}_{N(\mathbf{s}_i)}$ 的协方差矩阵，让 $\mathbf{C}_{\mathbf{s}_i, N(\mathbf{s}_i)}$ 为随机向量 $\mathbf{w}(\mathbf{s}_i)$ 和 $\mathbf{w}_{N(\mathbf{s}_i)}$ 之间大小为 $q×mq$ 的互协方差矩阵。标准分布理论揭示：

$$\tilde{p}(\mathbf{w}_{\mathscr{S}}) = \prod^{k}_{i=1} N \left (\mathbf{w}(\mathbf{s}_i) \mid \mathbf{B}_{\mathbf{s}_i} \mathbf{w}_{N(\mathbf{s}_i)}, \mathbf{F}_{\mathbf{s}_i} \right ) \tag{3}$$

其中 $\mathbf{B}_{\mathbf{s}_i} = \mathbf{C}_{\mathbf{s}_i,N(\mathbf{s}_i)} \mathbf{C}_{N(\mathbf{s}_i)}^{-1}$ 和 $\mathbf{F}_{\mathbf{s}_i} = \mathbf{C}(\mathbf{s}_i, \mathbf{s}_i) − \mathbf{C}_{\mathbf{s}_i,N(\mathbf{s}_i)} \mathbf{C}_{N(\mathbf{s}_i)}^{-1} \mathbf{C}_{N(\mathbf{s}_i), \mathbf{s}_i}$。附录 A2（可在线获取）显示 式 (3) 中的 $\tilde{p}(\mathbf{w}_\mathscr{S})$ 是具有协方差矩阵 $\tilde{\mathbf{C}}_{\mathscr{S}}$ 的多元高斯密度，这显然不同于 $\mathbf{C}_{\mathscr{S}}$。此外，如果 $N(\mathbf{s}_i)$ 对于 $\mathscr{S}$ 中每个 $\mathbf{s}_i$ 最多有 $m$ 个成员（ $m \ll k$ ），则 $\tilde{\mathbf{C}}^{−1}_{\mathscr{S}}$ 是稀疏的，最多有 $km(m+ 1)q^2/2$ 个非零元素。因此，对于某类非常普遍的邻居集，式 (2) 中定义的 $\tilde{p}(\mathbf{w}_{\mathscr{S}})$ 是具有稀疏精度矩阵的多元高斯分布的联合密度。

转向邻居集，选择 $N(\mathbf{s}_i)$ 为 $\{\mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_{i-1}\}$ 的任何子集保证了一个无环图 $\mathscr{G}$ 以及 式 (3) 中的一个有效概率密度。在 似然近似（likelihood approximation） 的上下文中存在几种特殊情况。例如，Vecchia (1988)^[46] 和 Stroud、Stein 和 Lysen (2014)^[45] 将 $N(\mathbf{s}_i)$ 指定为集合 $\{\mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_{i-1}\}$ 中的 $m$ 个 $\mathbf{s}_i$ 的欧式最近邻。Stein、Chi 和 Welty (2004) ^[44] 考虑了 $\{\mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_{i-1}\}$ 中的最近和最远邻居。Gramacy 和 Apley (2014)^[21] 在选择 $N(\mathbf{s}_i)$ 时提供了更大的灵活性，但可能需要多个近似值才能有效。

以上所有选择都取决于位置的排序。空间位置不是自然排序的，可以强加顺序，例如，对坐标之一进行排序。当然，可以使用坐标的任何其他函数来强加顺序。上述作者已经有力地证明了排序的选择对 式 (1) 与 式 (3) 的近似没有明显影响。我们自己的模拟实验（参见附录 A5，可在线获取）与这些发现一致；基于 $\tilde{p}(\mathbf{w}_\mathscr{S})$ 的推断对于位置排序而言非常稳健。这并不完全令人惊讶。显然，无论我们在 式 (1) 中选择什么顺序，$p(\mathbf{w}_{\mathscr{S}})$ 都会产生完整的联合密度。请注意，我们根据根据 式 (1) 中的特定顺序构造的邻居集将 式 (1) 简化为 式 (2)。 式 (1) 中的不同排序将为 式 (2) 产生一组不同的邻居。由于 $\tilde{p}(\mathbf{w}_\mathscr{S})$ 最终依赖于从邻居那里借来的信息，其有效性通常取决于邻居数量而不是特定顺序。

在下一节中，我们将把近似密度 $\tilde{p}(\mathbf{w}_{\mathscr{S}})$ 扩展到一个空间过程。我们的后续扩展适用于能够确保无环的任何 $N(\mathbf{s}_i)$ 选择。一般来说，确定能够最佳预测 $\mathbf{s}_i$ 的 $m$ 个位置构成的 “最佳子集” ，本身是一个非凸优化问题。这很难实现并且违背了使用较小条件集来简化计算的目的。尽管如此，我们发现 Vecchia 从 $\{\mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_{i-1}\}$ 中选择 $m$ 个最近邻的方法很简单，并且在广泛的模拟实验中表现得非常好。接下来，这也将是我们对 $N(\mathbf{s}_i)$ 的选择，相应的密度 $\tilde{p}(\mathbf{w}_{\mathscr{S}})$ 将被称为 $\mathbf{w}_{\mathscr{S}}$ 的 “最近邻密度”。

2.2 扩展到高斯过程

令 $\mathbf{u}$ 是位于 $\mathcal{D}$ 内 $\mathscr{S}$ 外的任何位置。与 $N(\mathbf{s}_i)$ 的定义一致，设 $N(\mathbf{u})$ 为 $\mathbf{u}$ 中的 $m$ 个最近邻的集合。则对于任何有限集 $\mathscr{U} = \{\mathbf{u}_1, \mathbf{u}_2, \ldots , \mathbf{u}_r\}$，$\mathscr{U} \cap \mathscr{S}$ 为空，我们定义以 $\mathbf{w}_{\mathscr{S}}$ 为条件的 $\mathbf{w}_{\mathscr{U}}$ 的最近邻密度为：

$$\tilde{p}(\mathbf{w}_{\mathscr{U}} \mid \mathbf{w}_{\mathscr{S}}) = \prod^{r}_{i=1} p \left(\mathbf{w}(\mathbf{u}_i) \mid \mathbf{w}_{N(\mathbf{u}_i)} \right) \tag{4}$$

这个条件密度类似于式 (2)（除了所有邻居集都是 $\mathscr{S}$ 的子集），这确保了正确的条件密度。实际上式 (2) 和式 (4) 足以描述域 $\mathcal{D}$ 上任意有限集的联合密度。更准确地说，如果 $\mathscr{V}= \{ \mathbf{v}_1, \mathbf{v}_2, \ldots , \mathbf{v}_n\}$ 是 $\mathcal{D}$ 中的任何有限子集，则使用式 (4) 我们获得 $\mathbf{w}_{\mathscr{V}}$ 的密度为:

$$\tilde{p}(\mathbf{w}_{\mathscr{V}}) = \int \tilde{p}(\mathbf{w}_{\mathscr{U}} \mid \mathbf{w}_{\mathscr{S}}) \tilde{p}(\mathbf{w}_{\mathscr{S}}) \prod_{\{\mathbf{s}_i \in \mathscr{S}-\mathscr{V} \}} d(\mathbf{w}(\mathbf{s}_i)) \text{ , where } \mathscr{U} = \mathscr{V} - \mathscr{S} \tag{5}$$

如果 $\mathscr{U}$ 为空，则式 (4) 意味着式 (5) 中的 $\tilde{p}(\mathbf{w}_{\mathscr{U}} \mid \mathbf{w}_{\mathscr{S}}) = 1$。如果 $\mathscr{S}-\mathscr{V}$ 为空，则不需要式 (5) 中的积分。

这些在有限拓扑上定义的概率密度符合 Kolmogorov 一致性准则，因此对应于一个 $\mathcal{D}$ 上的有效空间过程（参见附录 A3，可在线获取）。因此，给定任何原始（父）空间过程和任何固定参考集 $\mathscr{S}$ ，我们可以使用 $\mathscr{S}$ 中的邻居集集合在域 $\mathcal{D}$ 上构造一个新过程。我们将此过程称为从原始父过程派生的 “最近邻过程”。如果父过程是 $\mathcal{GP}(\mathbf{0}, \mathbf{C}(·,· \mid \boldsymbol{\theta}))$，则：

$$\tilde{p}(\mathbf{w}_{\mathscr{U}} \mid \mathbf{w}_{\mathscr{S}}) = \prod^{r}_{i=1} N(\mathbf{w}(\mathbf{u}_i) \mid \mathbf{B}_{\mathbf{u}_i} \mathbf{w}_{N(\mathbf{u}_i)}, \mathbf{F}_{\mathbf{u}_i}) = N(\mathbf{B}_{\mathscr{U}} \mathbf{w}_{\mathscr{S}}, \mathbf{F}_{\mathscr{U}} ) \tag{6}$$

对于 $\mathcal{D}$ 内 $\mathscr{S}$ 外的任何有限集 $\mathscr{U}= \{\mathbf{u}_1, \mathbf{u}_2, \ldots , \mathbf{u}_r\}$，其中 $\mathbf{B}_{\mathbf{u}_i}$ 和 $\mathbf{F}_{\mathbf{u}_i}$ 的定义类似于式 (3)，基于邻居集 $N(\mathbf{u}_i)$，$\mathbf{F}_{\mathscr{U}} = \text{diag}(F_{\mathbf{u}_1}, F_{\mathbf{u}_2}, \ldots , F_{\mathbf{u}_r})$，$\mathbf{B}_{\mathscr{U}}$ 是一个大小为 $nq×kq$ 的稀疏矩阵，每行最多有 $mq$ 个非零元素（参见附录 A4，可在线获取）。

对于 $\mathcal{D}$ 中的任何有限集 $\mathscr{V}$，$\tilde{p}(\mathbf{w}_{\mathscr{V}})$ 是$\mathscr{V}$ 上具有如下互协方差函数的高斯过程实现的密度。

$$\widetilde{\mathbf{C}}\left(\mathbf{v}_1, \mathbf{v}_2 ; \boldsymbol{\theta}\right)= \begin{cases} \widetilde{\mathbf{C}}_{\mathbf{s}_i, \mathbf{s}}, \quad \text { if } \mathbf{v}_1=\mathbf{s}_i \text { and } \mathbf{v}_2=\mathbf{s}_j \text { are both in } \mathcal{S}, \\ \mathbf{B}_{\mathbf{v}_1} \widetilde{\mathbf{C}}_{N\left(\mathbf{v}_1\right), \mathbf{s}_j} & \text { if } \mathbf{v}_1 \notin \mathcal{S} \text { and } \mathbf{v}_2=\mathbf{s}_j \in \mathcal{S}, \\ \mathbf{B}_{\mathbf{v}_1} \widetilde{\mathbf{C}}_{N\left(\mathbf{v}_1\right), N\left(\mathbf{v}_2\right)} \mathbf{B}_{\mathbf{v}_2}^{\prime}+\delta \delta_{\left(\mathbf{v}_1=\mathbf{v}_2\right)} \mathbf{F}_{\mathbf{v}_1}, \quad \text { if } \mathbf{v}_1 \text { and } \mathbf{v}_2 \text { are noth in } \mathcal{S} \end{cases}$$

其中 $\mathbf{v}_1$ 和 $\mathbf{v}_2$ 是 $\mathcal{D}$ 中的任意两个位置，$\tilde{\mathbf{C}}_{A,B}$ 表示由集合 $A$ 和 $B$ 中的位置索引的子矩阵 $\tilde{\mathbf{C}}_{\mathscr{S}}$ ，$δ_{(\mathbf{v}_1=\mathbf{v}_2)}$ 是 Kronecker delta。附录 A4（可在线获取）还表明 $\tilde{\mathbf{C}}(\mathbf{v}_1, \mathbf{v}_2 | \boldsymbol{\theta})$ 对于勒贝格测度零集合之外的所有 $(\mathbf{v}_1, \mathbf{v}_2)$ 对，都是连续的。

这就完成了从父高斯过程 $\mathcal{GP}(\mathbf{0}, \mathbf{C}(·, · | \boldsymbol{\theta}))$ 导出最近邻高斯过程 $\mathcal{NNGP}(\mathbf{0}, \tilde{\mathbf{C}}(·, · | \boldsymbol{\theta}))$ 的构造。在最近邻高斯过程中，$k$ 的大小可以与数据集大小一致，甚至更大。计算复杂度的降低是通过最近邻高斯过程的精度矩阵中的稀疏性实现的。与低秩过程不同，最近邻高斯过程不是退化过程。它是一个正确的、稀疏性归纳的高斯过程，可立即作为分层建模的先验使用，并且正如在下一节中展示的那样，提供了巨大的计算优势。

3. 贝叶斯估计与实现

3.1 分层模型

考虑一个 $l$ 个因变量的向量 $y(\mathbf{t})$，在位置 $\mathbf{t} \in \mathcal{D} \subseteq \mathbb{R}^d$ 处的一个空间变化回归模型：

$$y(\mathbf{t}) = \mathbf{X}(\mathbf{t})^{\prime}\boldsymbol{\beta} + \mathbf{Z}(\mathbf{t})^{\prime}\mathbf{w}(\mathbf{t}) + \boldsymbol{\varepsilon}(\mathbf{t}) \tag{8}$$

其中 $\mathbf{X}(\mathbf{t})^{\prime}$ 是大小为 $l× p$ 的固定空间位置预测变量矩阵，$\mathbf{w}(\mathbf{t})$ 是一个 $q × 1$ 空间过程，该过程形成大小为 $l × q$ 的固定设计矩阵 $\mathbf{Z}(\mathbf{t})^{\prime}$ 的系数，$\boldsymbol{\varepsilon}(\mathbf{t}) \stackrel{i.i.d.}{\sim} \mathcal{N}(\mathbf{0, D})$ 是一个 $l × 1$ 的白噪声过程，它捕获测量误差或具有扩散矩阵 $D$ 的微尺度可变性，我们假设它是对角矩阵，元素为 $τ_j^2, j= 1, 2, \ldots , l$。矩阵 $\mathbf{X}(\mathbf{t})^{\prime}$ 是 $p = \sum^l_{i =1} p_i$ 的块对角矩阵，其中 $1 × p_i$ 的向量 $\mathbf{x}_i(\mathbf{t})^{\prime}$ （可能包括截距）是第 $i$ 个块（对于每一个 $i= 1, 2,\ldots, l$）。式 (8) 中的模型蕴含了几个特定的空间模型。例如，让 $q = l$ 和 $\mathbf{Z}(\mathbf{t})^{\prime} = \mathbf{I}_{l×l}$ 会形成一个多变量空间回归模型，其中 $\mathbf{w}(\mathbf{t})$ 充当空间变化的截距。另一方面，我们可以设想所有系数都在随空间变化，并设置 $q= p$ 且 $\mathbf{Z}(\mathbf{t})^{\prime} = \mathbf{X}(\mathbf{t})^{\prime}$。

对于可扩展性来说，可以替换掉式 (8) 中 $\mathbf{w}(\mathbf{t})$ 的常规高斯过程先验，我们假设 $\mathbf{w}(\mathbf{t}) \sim \mathcal{NNGP}(0, \tilde{\mathbf{C}}(·, · | \boldsymbol{\theta}))$ 是从父过程 $\mathcal{GP}(0, \mathbf{C}( ·, · | \boldsymbol{\theta}))$ 派生出来的。任何有效的各向同性互协方差函数都可用于构造 $\mathbf{C}(·, · | \boldsymbol{\theta})$（参见例如 Gelfand 和 Banerjee 2010 ^[17]）。为了阐明，设 $\mathscr{T}= \{\mathbf{t}_1, \mathbf{t}_2, \ldots , \mathbf{t}_n\}$ 是已观测响应变量和预测变量的位置集。该集合可以但不必与最近邻高斯过程的参考集合 $\mathscr{S}= \{\mathbf{s}_1, \mathbf{s}_2, \ldots , \mathbf{s}_k\}$ 相交。不失一般性，我们将 $\mathscr{T}$ 分成 $\mathscr{S}^*$ 和 $\mathscr{U}$ ，其中 $\mathscr{S}^{*} = \mathscr{S} \cap \mathscr{T} = \{\mathbf{s}_{i1}, \mathbf{s}_{i2}, \ldots , \mathbf{s}_{ir}\}$，$\mathbf{s}_{i_j}= \mathbf{t}_j,\quad j= 1, 2, \ldots , r$ 且 $\mathscr{U} =\mathscr{T} - \mathscr{S} = \{ \mathbf{t}_{r+1} , \mathbf{t}_{r+2}, \ldots, \mathbf{t}_n\}$。由于 $\mathscr{S} \cup \mathscr{T} = \mathscr{S} \cup \mathscr{U}$，我们可以根据父过程在 $\mathscr{S}$ 和 $\mathscr{U}$ 上的实现来完全指定最近邻高斯过程的实现，并分层地表示为 $\mathbf{w}_{\mathscr{U}} \mid \mathbf{w}_{\mathscr{S}} \sim \mathcal{N}(\mathbf{B}_{\mathscr{U}} \mathbf{w}_{\mathscr{S}} , \mathbf{F}_{\mathscr{U}})$ 和 $\mathbf{w}_{\mathscr{S}} \sim \mathcal{N}(\mathbf{0}, \tilde{\mathbf{C}}_{\mathscr{S}})$。对于完整的贝叶斯规格，我们进一步指定了 $\boldsymbol{\beta}$、$\boldsymbol{\theta}$ 和 $τ_j^2$ 的先验分布。例如，根据惯例先验指定，我们可以得到联合分布：

$$p(\boldsymbol{\theta}) × \prod^{l}_{j=1} \mathcal{IG}(\tau^{2}_{j} \mid a_{\tau_j}, b_{\tau_j}) × N(\boldsymbol{\beta} \mid \boldsymbol{μ}_{\boldsymbol{\beta}}, \mathbf{V}_{\boldsymbol{\beta}}) × N(\mathbf{w}_{\mathscr{U}} \mid \mathbf{B}_{\mathscr{U}} \mathbf{w}_{\mathscr{S}}, \mathbf{F}_{\mathscr{U}}) × N(\mathbf{w}_{\mathscr{S}} \mid \mathbf{0}, \tilde{\mathbf{C}}_{\mathscr{S}}) × \prod^{n}_{i=1} N(\mathbf{y}(\mathbf{t}_i) \mid \mathbf{X}(\mathbf{t}_i)^{\prime} \boldsymbol{\beta} + \mathbf{Z}(\mathbf{t}_i)^{\prime} \mathbf{w}(\mathbf{t}_i), \mathbf{D}) \tag{9}$$

其中 $p(\boldsymbol{\theta})$ 是 $\boldsymbol{\theta}$ 上的先验，$\mathcal{IG}(\tau^{2}_{j} \mid a_{\tau_j}, b_{\tau_j})$ 表示逆 Gamma 密度。

3.2 参数估计与空间预测

为了描述用于估计 式 (9) 的 Gibbs 采样器，我们类似地定义 $\mathbf{y} = (\mathbf{y}(\mathbf{t}_1)^{\prime}, \mathbf{y}(\mathbf{t}_2)^{\prime}, \ldots , \mathbf{y}(\mathbf{t}_n)^{\prime})^{\prime}$，以及 $\mathbf{w}$ 和 $\boldsymbol{\varepsilon}$。此外，引入 $\mathbf{X} = [\mathbf{X}(\mathbf{t}_1) : \mathbf{X}(\mathbf{t}_2) : \ldots : \mathbf{X}(\mathbf{t}_n)]'$, $\mathbf{Z} = \text{diag}(\mathbf{Z}(\mathbf{t}_1)^{\prime}, \ldots , \mathbf{Z}(\mathbf{t}_n)^{\prime})$，以及 $\mathbf{D}_n= \text{Cov}(\boldsymbol{\varepsilon}) = \text{diag}(D, \ldots, D)$。

（1）固定效应参数的估计

$\boldsymbol{\beta}$ 的完整条件分布是 $\mathcal{N}(\mathbf{V}_{\boldsymbol{\beta}}^{∗} \boldsymbol{μ}_{\boldsymbol{\beta}}^{∗}, \mathbf{V}_{\boldsymbol{\beta}}^{∗})$，其中 $\mathbf{V}_{\boldsymbol{\beta}}^{∗} = (\mathbf{V}_{\boldsymbol{\beta}}^{-1} + \mathbf{X}^{\prime} \mathbf{D}_{n}^{−1} \mathbf{X})^{−1}$, $\boldsymbol{μ}_{\boldsymbol{\beta}}^{∗ }= (\mathbf{V}_{\boldsymbol{\beta}}^{−1} \boldsymbol{μ}_{\boldsymbol{\beta}} + \mathbf{X}^{\prime} \mathbf{D}_{n}^{−1} (\mathbf{y − Zw}))$。

（2）噪声超参数的估计

$\tau^{2}_{j}$ 的逆伽马先验导致共轭完全条件分布 $\mathcal{IG}(a_{\tau_j} + \frac{n}{2}, b_{\tau_j} + \frac{1}{2} (\mathbf{y}_{∗ j} − \mathbf{X}_{∗j} \boldsymbol{\beta} − \mathbf{Z}_{∗j} \mathbf{w})^{\prime} (\mathbf{y}_{*j} − \mathbf{X}_{*j} \boldsymbol{\beta} − \mathbf{Z}_{*j} \mathbf{w})$，其中 $y^{*j}$ 是指包含 $y(\mathbf{t}_i)$ 的第 $j$ 个坐标的 $n×1$ 向量，$\mathbf{X}^{*j}$ 和 $Z^{*j}$ 分别是相应的固定和空间效应协变量矩阵。

（3）协方差超参数的估计

为了更新 $\boldsymbol{\theta}$，我们使用目标密度为 $p(\boldsymbol{\theta}) × N(\mathbf{w}_{\mathscr{S}} \mid \mathbf{0}, \tilde{\mathbf{C}}_{\mathscr{S}}) × N(\mathbf{w}_{\mathscr{U}} \mid \mathbf{B}_{\mathscr{U}} \mathbf{w}_{\mathscr{S}} , \mathbf{F}_{\mathscr{U}})$ 的随机游走 Metropolis 步骤，其中:

$$\begin{align*} &N(\mathbf{w}_{\mathscr{S}} \mid \mathbf{0}, \tilde{\mathbf{C}}_{\mathscr{S}}) = \prod^{k}_{i=1} N(\mathbf{w}(\mathbf{s}_i) \mid \mathbf{B}_{\mathbf{s}_i} \mathbf{w}_{N(\mathbf{s}_i)}, \mathbf{F}_{\mathbf{s}_i})\\ &N(\mathbf{w}_{\mathscr{U}} \mid \mathbf{B}_{\mathscr{U}} \mathbf{w}_{\mathscr{S}}, \mathbf{F}_{\mathscr{U}}) = \prod^{n}_{i=r+1} N(\mathbf{w} (\mathbf{t}_i) \mid \mathbf{B}_{\mathbf{t}_i} \mathbf{w}_{N(\mathbf{t}_i)}, \mathbf{F}_{\mathbf{t}_i}) \end{align*} \tag{10}$$

可以在不进行任何 $n$ 维矩阵运算的情况下评估式 (10) 右侧乘积符号下的每个分量密度。

由于 $\mathbf{w}_{\mathscr{U}} \mid \mathbf{w}_{\mathscr{S}}$ 的组分是独立的，我们可以从其完整条件 $N(\mathbf{V}_{\mathbf{t}_i} \boldsymbol{μ}_{\mathbf{t}_i}, \mathbf{V}_{\mathbf{t}_i})$ 更新 $\mathbf{w}(\mathbf{t}_i)$，其中 $i= r+ 1, r+ 2,\ldots, n$， $\mathbf{V}_{\mathbf{t}_i} = (\mathbf{Z}(\mathbf{t}_i) \mathbf{D}^{−1} \mathbf{Z}(\mathbf{t}_i)^{\prime} + \mathbf{F}_{\mathbf{t}_i}^{−1})^{−1}$ 和 $\boldsymbol{μ}_{\mathbf{t}_i} = \mathbf{Z}(\mathbf{t}_i) \mathbf{D}^{−1} (\mathbf{y}(\mathbf{t}_i) − \mathbf{X}(\mathbf{t}_i)^{\prime} \boldsymbol{\beta}) + \mathbf{F}_{\mathbf{t}_i}^{−1} \mathbf{B}_{\mathbf{t}_i} \mathbf{w}_{N(\mathbf{t}_i)}$。最后，我们逐一更新 $\mathbf{w}_{\mathscr{S}}$ 的组件。

对于 $\mathcal{D}$ 中的任意两个位置 $\mathbf{s}$ 和 $\mathbf{t}$，如果 $\mathbf{s} \in N(\mathbf{t})$ 并且是 $N(\mathbf{t})$ 的第 $l$ 个分量，即 $\mathbf{s} = N(\mathbf{t})(l)$，则将 $\mathbf{B_{t,s}}$ 定义为 $l× l$ 的子矩阵，由 $\mathbf{B_t}$ 的列 $(l−1)q+ 1, (l−1)q + 2, \ldots , lq$ 组成。让 $U(\mathbf{s}_i) = \{\mathbf{t} \in \mathscr{S} \cup \mathscr{T} \mid \mathbf{s}_i \in N(\mathbf{t})\}$ ，并且对于每个 $\mathbf{t} \in U(\mathbf{s}_i)$ 定义 $\mathbf{a}_{\mathbf{t},\mathbf{s}_i} =\mathbf{w}(\mathbf{t}) − \Sigma_{s \in N(\mathbf{t}),\mathbf{s} \neq \mathbf{s}_{i}} \mathbf{B}_{\mathbf{t,s}} \mathbf{w}(\mathbf{s})$。然后，对于 $i= 1, 2, \ldots , k$，我们有完整的条件 $\mathbf{s}_{\mathbf{s}_i} \mid · \sim \mathcal{N}(\mathbf{V}_{\mathbf{s}_i} \boldsymbol{\mu}_{\mathbf{s}+i}, \mathbf{V}_{\mathbf{s}_i})$，其中

$$\begin{align*} &\mathbf{V}_{\mathbf{s}_i} = (\text{In}(\mathbf{s}_i \in \mathscr{S}^∗) \mathbf{Z}(\mathbf{s}_i) \mathbf{D}^{−1} \mathbf{Z}(\mathbf{s}_i)^{\prime} + \mathbf{F}_{\mathbf{s}_i}^{ −1} + \Sigma_{\mathbf{t} \in U(\mathbf{s}_i)} \mathbf{B}_{\mathbf{t}, \mathbf{s}_i}^{\prime} \mathbf{F}_{\mathbf{t}}^{−1} \mathbf{B}_{\mathbf{t}, \mathbf{s}_i})^{−1} \\ &\boldsymbol{\mu}_{\mathbf{s}_i} = \text{In}(\mathbf{s}_i \in \mathscr{S}^∗) \mathbf{Z}(\mathbf{s}_i) \mathbf{D}^{−1} (\mathbf{y}(\mathbf{s}_i) − \mathbf{X}(\mathbf{s}_i)^{\prime} \boldsymbol{\beta}) + \mathbf{F}_{\mathbf{s}_i}^{−1} \mathbf{B}_{\mathbf{s}_i} \mathbf{w}_{N(\mathbf{s}_i)} + \Sigma_{\mathbf{t} \in U (\mathbf{s}_i)} \mathbf{B}_{\mathbf{t}, \mathbf{s}_i}^{\prime} \mathbf{F}_{\mathbf{t}}^{−1} \mathbf{a}_{\mathbf{t},\mathbf{s}_i} \end{align*}$$

其中 $\text{In}(·)$ 表示指示函数。因此，$w^{\prime}$ 也可以在不需要存储或分解任何 $n × n$ 矩阵的情况下更新。

（4）空间预测

假设 $\mathbf{t}$ 是我们打算在给定 $\mathbf{X}(\mathbf{t})$ 和 $\mathbf{Z}(\mathbf{t})$ 的情况下预测 $\mathbf{y}(\mathbf{t})$ 的新位置。用于估计的 Gibbs 采样器还生成后验样本 $\mathbf{w}_{\mathscr{S} \cup \mathscr{T}} \mid \mathbf{y}$。所以，如果 $\mathbf{t} \in \mathscr{S} \cup \mathscr{T}$ ，那么我们简单地从 $N(\mathbf{X}(\mathbf{t})^{\prime}\boldsymbol{\beta} + \mathbf{Z}(\mathbf{t})^{\prime}\mathbf{w}(\mathbf{t}), \mathbf{D})$ 得到 $\mathbf{y}(\mathbf{t})|\mathbf{y}$ 的样本。如果 $t$ 在 $\mathscr{S} \cup \mathscr{T}$ 之外，那么我们从其完整条件 $N(\mathbf{B_t} \mathbf{w}_{\mathscr{S}}, \mathbf{F_t})$ 生成 $\mathbf{w}(\mathbf{t})$ 的样本，随后生成 $\mathbf{y}(\mathbf{t}) \mid \mathbf{y}$ 的后验样本。

3.3 计算复杂度

实施 第 3.2 节 中的最近邻高斯过程模型表明，无需任何大型矩阵运算即可完成 Gibbs 采样器的整个传递。式 (9) 和完整的地统计分层模型之间的唯一区别是空间过程被建模为最近邻高斯过程先验而不是标准高斯过程。为了进行比较，我们提供了 flops 计数的粗略估计，以在采样器的每次迭代中生成 $\boldsymbol{\theta}$ 和 $\mathbf{w}$。我们仅根据样本大小 $n$、参考集大小 $k$ 和邻居集大小 $m$ 来表示计算复杂度，假设其他维度很小。对于所有位置，$\mathbf{t} \in \mathscr{S} \cup \mathscr{T}$ 、$\mathbf{B_t}$ 和 $\mathbf{F_t}$ 可以使用 $\mathcal{O}(m^3)$ flops 来计算。因此，从式 (10) 中很容易看出可以使用 $\mathcal{O}((n+k)m^3)$ 次 flops 计算 $p(\boldsymbol{\theta}|·)$。为 $\mathbf{w}$ 和 $\boldsymbol{\theta}$ 生成一组后验样本的所有后续计算都需要大约 $\mathcal{O}((n+k) m^2)$ 次 flops 。

因此，总 flops 的阶数为 $(n+k)m^3$，因此与 $\mathscr{S} \cup \mathscr{T}$ 中的位置总数成线性关系。这确保了最近邻高斯过程对大型数据集的可扩展性。将其与具有密集相关矩阵的完整高斯过程模型进行比较，后者需要 $\mathcal{O}(n^3)$ flops 来在每次迭代中更新 $\mathbf{w}$。第 5.1 节 和在线附录 A6 中的模拟结果表明，通常具有非常小的 $m(\approx 10)$ 值的最近邻高斯过程模型提供的推断几乎与完整的地统计模型无法区分。因此，对于较大的 $n$，flops 计数会大大减少，并且相对于 $k$ 是线性的，从而保证了即使 $k \approx n$ 也存在可行实现。

这提供了相对于低秩模型的可扩展性，在低秩模型中，计算成本是 “结” 数量的二次方，从而限制了结集的大小。此外，完整的地统计模型需要存储 $n × n$ 的距离矩阵，这可能会耗尽大型数据集的存储资源。最近邻高斯过程模型只需要每个位置的邻居之间的距离矩阵，从而存储 $n+ k$ 小矩阵，每个矩阵的阶数为 $m× m$。

3.4 模型比较与超参数的选择

如第 2 节所述，给定任何父高斯过程和任何固定参考位置集，我们可以构建有效的最近邻高斯过程。最近邻高斯过程的有限维似然最终取决于参考集 $\mathscr{S}$ 的选择和每个 $N(\mathbf{s}_i)$ 的大小，即 $m$。

（1）参考点集 $\mathcal{S}$ 的选择

选择参考集类似于为预测过程选择结点。与低秩模型中的“结”数量不同，$\mathscr{S}$ 的大小不会阻碍计算的可扩展性。由于最近邻高斯过程模型中的 flop 计数仅随 $\mathscr{S}$ 的大小线性增加，因此 $\mathscr{S}$ 中的位置数量可以很大，对于 $\mathscr{S}$ 的选择可以更灵活。

跨越整个域的网格上的点似乎是 $\mathscr{S}$ 的合理选择。例如，我们可以使用密集网格构建一个大型 $\mathscr{S}$ 来提高性能，而不会对计算成本产生不利影响。对于大型数据集，另一个可能更简单的选项是简单地固定 $\mathscr{S} =\mathscr{T}$ ，即观测到的位置集。这种选择通过避免在 Gibbs 采样器中对 $\mathbf{w}$ 进行额外采样，进一步降低了计算成本。我们的实证调查（参见 第 5.1 节）表明，选择 $\mathscr{S}=\mathscr{T}$ 提供的推断与选择 $\mathscr{S}$ 为大型数据集域上的网格几乎没有区别。

Stein、Chi 和 Welty (2004)^[44] 以及 Eidsvik 等 (2014) ^[12] 提出使用三明治方差估计器来评估基于邻居的伪似然推断能力。 Shaby (2012) ^[36] 开发了一种采样后的三明治方差调整，用于使用伪似然法的准贝叶斯方法参数的后验可信区间。然而，用于获得三明治方差估计量的渐近结果是基于难以在具有不规则放置的数据点的空间设置中验证的假设。

此外，我们将最近邻高斯过程视为用于拟合数据的独立模型，而不是原始高斯过程的近似值。因此，我们避免进行这种三明治方差调整。相反，我们可以简单地使用任何标准模型比较指标，例如偏差信息标准（DIC；Spiegelhalter 等，2002 年 ^[38]）、GPD 分数（Gelfand 和 Ghosh，1998 年 ^[19]）或均方根预测误差/变异系数的均方根误差（RMSPE/RMSECV；Yeniay 和 Goktas 2002 ^[50]）比较最近邻高斯过程和任何其他候选模型的性能。

（2）邻居数量 $m$ 的选择

相同的模型比较指标也用于选择 $m$。然而，正如我们稍后在第 5.1 节中说明的那样，通常 $10$ 到 $15$ 之间的小值 $m$ 会产生与完整地统计模型相当的性能。虽然较大的 $m$ 可能有利于海量数据集，但也许在不同设计方案下，它会比可比较推断的低秩模型中的结数小得多（参见第 5.1 节）。

4 改进的最近邻高斯过程模型和算法

4.1 使用稀疏 Cholesky 对 w 进行块更新

第 3.2 节 中详述的 Gibbs 采样算法对于每次迭代具有线性 flop 计数的大型数据集非常有效。但是，由于 $\mathbf{w}_{\mathscr{S}}$ 中元素的顺序更新，它有时会遇到收敛速度慢的问题。顺序更新的替代品是对 $\mathbf{w}_{\mathscr{S}}$ 进行块更新。我们选择 $\mathscr{S} =\mathscr{T}$ 使得对于所有 $i= 1, 2, \ldots , n$，$\mathbf{s}_i = \mathbf{t}_i$；并且我们用 $\mathbf{w}$ 表示 $\mathbf{w}_{\mathscr{S}} = \mathbf{w}_{\mathscr{T}}$。然后，

$$\mathbf{w} \mid \cdot \sim \mathcal{N}(\mathbf{V}_{\mathscr{S}} \mathbf{Z}^{\prime} \mathbf{D}_n^{-1} (\mathbf{y} − \mathbf{X} \boldsymbol{\beta}), \mathbf{V}_{\mathscr{S}}) \tag{11}$$

其中 $\mathbf{V}_{\mathscr{S}} = (\mathbf{Z}^{\prime} \mathbf{D}_n^{-1} \mathbf{Z} + \tilde{\mathbf{C}}_{\mathscr{S}}^{ −1})^{−1}$ 。

回想一下 $\tilde{\mathbf{C}}_{\mathscr{S}}^{-1}$ 是稀疏的。由于 $\mathbf{Z}$ 和 $\mathbf{D}_n$ 是块对角的，$\mathbf{V}^{-1}_{\mathscr{S}}$ 保留了 $\tilde{\mathbf{C}}^{-1}_{\mathscr{S}}$ 的稀疏性。因此， $\tilde{\mathbf{C}}^{-1}_{\mathscr{S}}$ 的稀疏 Cholesky 分解将有效地产生 $\mathbf{V}_{\mathscr{S}}$ 的 Cholesky 因子。这将促进 Gibbs 采样器中 $\mathbf{w}$ 的块更新。

4.2 响应的最近邻高斯过程模型

另一种可能的方法涉及响应 $\mathbf{y} (\mathbf{s})$ 的最近邻高斯过程模型。

如果 $\mathbf{w}(\mathbf{s})$ 是高斯过程，那么 $\mathbf{y}(\mathbf{s}) = \mathbf{Z}(\mathbf{s})^{\prime}\mathbf{w}(\mathbf{s}) + \boldsymbol{\varepsilon}$ 也是（不失一般性，我们假设 $\boldsymbol{\beta}= 0$）。可以直接将最近邻高斯过程用于 $\mathbf{y}(\mathbf{s})$ 而不是 $\mathbf{w}(\mathbf{s})$。也就是说，我们从父高斯过程 $\mathcal{GP}(\mathbf{0}, \boldsymbol{\Sigma}(·, · | \boldsymbol{\theta}))$ 导出 $\mathbf{y}(\mathbf{s}) \sim \mathcal{NNGP}(\mathbf{0}, \tilde{\boldsymbol{\Sigma}}(·, ·))$。类似于第 3 节的 Gibbs 采样器现在享有避免 $\mathbf{w}$ 的完整条件的额外优势。这导致 Vecchia (1988)^[46] 和 Stein、Chi 和 Welty (2004)^[44] 的贝叶斯模拟，但排除了对空间残差表面 $\mathbf{w}(\mathbf{s})$ 的推断。对 $\mathbf{w}(\mathbf{s})$ 建模提供了对剩余空间表面的额外洞察，并且通常对于识别潜在协变量或引出无法解释的空间模式很重要。 Vecchia (1992) ^[47] 在空间模型上使用了最近邻近似值来观测除了通常的空间分量 $(\mathbf{w})$ 之外还具有独立测量误差 (nuggets) 的观测值 $(\mathbf{y})$。但是，使用这种方法可能无法恢复 $\mathbf{w}$。例如，对于某些 $\mathbf{w}(\mathbf{s}) \sim \mathcal{GP}(\mathbf{0}, \mathbf{C}(·,· | \boldsymbol{\theta}))$ 和白噪声过程 $\boldsymbol{\varepsilon}(\mathbf{s})$，具有块金效应的单变量平稳过程 $\mathbf{y(s)}$ 可以分解为 $\mathbf{y(s)} = \mathbf{w(s)} + \boldsymbol{ε}(\mathbf{s})$（令 $\boldsymbol{β}= \mathbf{0}$）。 如果 $\mathbf{y} = \mathbf{w} + \boldsymbol{\varepsilon}$，其中 $\mathbf{w} \sim N(\mathbf{0, C})$，$\boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \boldsymbol{τ}^2 \mathbf{I}_n)$，则 $\text{cov}(\mathbf{y}) = \mathbf{C} + \boldsymbol{τ}^2 \mathbf{I} = \boldsymbol{\Sigma}$，$\boldsymbol{\Sigma}$ 的所有特征值都大于 $\boldsymbol{τ}^2$，并且 $\text{cov}( \mathbf{w | y}) = \boldsymbol{\tau}^2 \mathbf{I}_n − \boldsymbol{τ}^4 \boldsymbol{\Sigma}^{−1}$。然而，对于 $\mathbf{y}(\mathbf{s}) \sim \mathcal{NNGP}(\mathbf{0}, \tilde{\boldsymbol{\Sigma}}(·, ·))$，$\tilde{\boldsymbol{\Sigma}}$ 的特征值可能小于 $\boldsymbol{\tau}^2$，因此对于每个 $\boldsymbol{\tau}^2 > 0$ 和 $p(\mathbf{w | y})$ 不再明确定义。

通过使用 $\mathbf{w}$ 的最近邻高斯过程先验，如式 (9) ，然后积分出 $\mathbf{w}$，可以获得一个不同的模型。由此产生的似然是 $\mathcal{N}(\mathbf{y} \mid \mathbf{X} \boldsymbol{\beta}, \boldsymbol{\Sigma}_{\mathbf{y}})$，其中 $\boldsymbol{\Sigma}_{\mathbf{y}}= \mathbf{Z} \tilde{\mathbf{C}}_{\mathscr{S}} \mathbf{Z}^{\prime} + \mathbf{D}_n$ 并且贝叶斯规范是使用 $\boldsymbol{\beta}$、$\tau^{2}_{j}$ 和 $\boldsymbol{\theta}$ 的先验完成的，如式 (9) 所示。该模型极大地减少了 Gibbs 采样器中的变量数量，同时保留了父模型中的块金效应。我们可以为边缘化模型中的参数生成如下所示的完整条件，：

$$\boldsymbol{\beta} \mid \mathbf{y},\boldsymbol{\phi} \sim \mathcal{N}(( \mathbf{V}_{\boldsymbol{\beta}}^{−1} + \mathbf{X}^{\prime} \boldsymbol{\Sigma}_\mathbf{y}^{-1} \mathbf{X})^{-1} (\mathbf{V}_{\boldsymbol{\beta}}^{-1} \boldsymbol{μ}_{\boldsymbol{\beta}} + \mathbf{X}^{\prime} \boldsymbol{\Sigma}_{\mathbf{y}}^{-1} \mathbf{y}), (\mathbf{V}_{\boldsymbol{\beta}}^{-1} + \mathbf{X}^{\prime} \boldsymbol{\Sigma}_{\mathbf{y}}^{-1} \mathbf{X})^{-1})。$$

很难从 $\boldsymbol{\Sigma}_{\mathbf{y}}^{-1}$ 中分解出 $\tau^{2}_{j}$ ，因此关于任何标准先验的共轭性都会丢失。对于任何易于处理的先验 $p(\boldsymbol{\theta})$，针对 $\boldsymbol{\theta}$ 的 Metropolis 块更新都是可行的。这涉及计算 $\mathbf{X}^{\prime} \boldsymbol{\Sigma}_{\mathbf{y}}^{-1}\mathbf{X}$、$\mathbf{X}^{\prime}\boldsymbol{\Sigma}_{\mathbf{y}}^{-1} \mathbf{y}$ 和 $(\mathbf{y} − \mathbf{X} \boldsymbol{\beta})^{\prime}\boldsymbol{\Sigma}_{\mathbf{y}}^{-1} (\mathbf{y} − \mathbf{X} \boldsymbol{\beta})$。

由于 $\boldsymbol{\Sigma}_{\mathbf{y}}^{-1} = \mathbf{D}_n^{-1} − \mathbf{D}_n^{-1} \mathbf{Z}(\tilde{\mathbf{C}}^{-1}_{\mathscr{S}} + \mathbf{Z}^{\prime} \mathbf{D}_n^{-1} \mathbf{Z})^{-1} \mathbf{Z}^{\prime} \mathbf{D}_n^{-1} = \mathbf{D}_n^{-1} − \mathbf{D}_n^{-1} \mathbf{Z} \mathbf{V}_{\mathscr{S}} \mathbf{Z}^{\prime} \mathbf{D}_n^{-1}$ ，其中 $\mathbf{V}_{\mathscr{S}}$ 由式 (11)给出，$\mathbf{V}^{−1}_{\mathscr{S}}$ 的稀疏 Cholesky 分解将是有益的。我们从 $p(\mathbf{w} \mid \mathbf{y}) = \int p(\mathbf{w} \mid \boldsymbol{\theta}, \boldsymbol{\beta}, \{\tau^{2}_{j}\}, \mathbf{y}) p(\boldsymbol{\theta}, \boldsymbol{\beta}, \{\tau^{2}_{j}\} \mid \mathbf{y})$ 使用组合采样为 $\mathbf{w}$ 抽取后验样本，我们从 $p(\mathbf{w} \mid \boldsymbol{\theta}^{(g)}, \boldsymbol{\beta}^{(g)}, \{\tau_{j}^{2(g)} \}, \mathbf{y})$ 中为每个采样的参数一对一抽取 $\mathbf{w}^{(g)}$。

在式 (9) 中对 $\mathbf{w}_{\mathscr{S}}$ 使用块更新和拟合 式 (9) 的边缘化版本，都需要一个有效的稀疏 Cholesky 求解器来求解 $V^{−1}_{\mathscr{S}}$。请注意，大多数稀疏 Cholesky 算法的计算费用取决于 $\tilde{\mathbf{C}}^{-1}_{\mathscr{S}}$ 的稀疏结构（主要是带宽）的精确性质（参见，例如，Davis 2006 ^[9]）。在这个边缘化模型中，吉布斯采样和预测所需的 flops 数量取决于 $\tilde{\mathbf{C}}^{-1}_{\mathscr{S}}$ 的稀疏结构，有时可能会大大超过非边缘化模型对 $\mathbf{w}_i$ 进行单独更新所实现的线性使用。因此，精确拟合算法的谨慎选择应基于给定数据集$\tilde{\mathbf{C}}^{-1}_{\mathscr{S}}$ 的稀疏结构。

4.3 时空和广义线性模型版本

在我们寻求离散时间点（例如，每周、每月或每年数据）的空间插值的时空场景中，我们将响应（可能是矢量值）写为 $\mathbf{y}_t(\mathbf{s})$，将随机效应写为 $\mathbf{w}_t(\mathbf{s})$。所需的推断包括每个时间点的空间插值。包含最近邻高斯过程的空间动态模型很容易表述如下：

$$\begin{align*} &\mathbf{y}_t (\mathbf{s}) = \mathbf{X}_t(\mathbf{s})^{\prime} \boldsymbol{\beta}_t + \boldsymbol{u}_t(\mathbf{s}) + \boldsymbol{\varepsilon}_t(\mathbf{s})\\ &\boldsymbol{\varepsilon}t(\mathbf{s}) \stackrel{i.i.d.}{\sim} \mathcal{N}(\mathbf{0}, \mathbf{D})\\ &\boldsymbol{\beta}_t = \boldsymbol{\beta}_{t − 1} + \boldsymbol{η}_t\\ &\boldsymbol{η}_t \stackrel{i.i.d.}{\sim} \mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma}_η)\\ &\boldsymbol{\beta}_0 \sim \mathcal{N}(\mathbf{m_0}, \boldsymbol{\Sigma_0})\\ &\mathbf{u}_t(\mathbf{s}) = \mathbf{u}_{t − 1} (\mathbf{s}) + \mathbf{w}_t(\mathbf{s})\\ &\mathbf{w}_t(\mathbf{s}) \stackrel{ind}{\sim} \mathcal{NNGP}(\mathbf{0}, \tilde{\mathbf{C}}( · , · \mid \boldsymbol{\theta}_t)) \end{align*}$$

因此，人们保留了与基于过程的空间动力学模型完全相同的结构，例如，与 Gelfand、Banerjee 和 Gamerman (2005) ^[18] 中的结构相同，并简单地将 $\mathbf{w}_t(\mathbf{s})$ 的独立高斯过程先验替换为独立的最近邻高斯过程，以实现计算易处理性。

以上说明了最近邻高斯过程在模型构建方面的吸引力和便利性。只需写下父模型，然后用最近邻高斯过程替换完整的高斯过程。作为一个定义明确的过程，最近邻高斯过程确保了有效的空间动态模型。类似地，可以构建动态时空卡尔曼滤波的最近邻高斯过程版本（例如，在 Wikle 和 Cressie 1999 中）^[49]。

使用空间广义线性模型处理非高斯（例如二值或计数）数据也很简单（GLM；Diggle、Tawn 和 Moyeed 1998 ^[10]；Lin 等人 2000 ^[29]；Kammann 和 Wand 2003 ^[25]；Banerjee、Carlin 和 Gelfand 2014 ^[1]）。在这里，最近邻高斯过程在第二阶段为随机效应提供了结构化依赖。

首先，我们用 $g(E(\mathbf{y}(\mathbf{t})))$ 替换 式 (8) 中的 $E[\mathbf{y}(\mathbf{t})]$ ，其中 $g(·)$ 是一个合适的链接函数，使得 $\boldsymbol{η}(\mathbf{t}) = g(E(\mathbf{y}(\mathbf{t}))) = \mathbf{X}(\mathbf{t})^{\prime} \boldsymbol{\beta} + \mathbf{Z}(\mathbf{t})^{\prime} \mathbf{w}(\mathbf{t})$。在第二阶段，我们将 $\mathbf{w}(\mathbf{t})$ 建模为最近邻高斯过程。第 3.2 节 和 第 3.3 节 中算法的优势仍然存在，但 第 4 节 中的某些替代算法可能不适用。例如，我们确实通过积分空间效应来获得易于处理的边缘化似然。

5 展示

参见原文。

6 总结与结论

我们将最近邻高斯过程视为大型地统计数据集的高度可扩展模型，而不是似然近似。

在推断性能和可扩展性方面，它明显优于高斯预测过程模型（Banerjee 等,2008 )^[3] 等竞争性低秩过程。参考集和生成的邻居集（大小为 $m$ ）定义了最近邻高斯过程。较大的 $m$ 会增加成本，但对于大型数据集，增加 $m$ 并没有明显的好处（参见 附录 A6，可在线获取）。虽然要求选择点的一些敏感性是预期的，但我们的结果表明推断非常稳健，并且 $m(< 20)$ 的非常适中的值通常就足够了。更大的数据集可能需要更大的参考集，但其大小不会妨碍计算。事实上，观测位置集是参考集最便利的选择。

这种选择的一个潜在问题是： 如果观测位置之间有很大的间隙（Gap），那么生成的最近邻高斯过程可能不是完整高斯过程的近似。这是因为参考集之外的观测会通过它们各自的邻居集相关，并且大的间隙可能意味着两个非常接近的点具有非常不同的邻居集，从而导致低相关性。我们在 附录 A7（可在线获取）中的模拟结果确实表明，在这种情况下，最近邻高斯过程的协方差场在间隙中的点处非常平坦。然而，即使选择了最近邻高斯过程模型，其性能也与完整的高斯过程模型相当，因为后者也无法提供关于较大间隙中观测的有力信息。当然，总是可以在整个域上选择一种网格作为 $\mathscr{S}$，来构造一个具有与完整高斯过程类似协方差函数的最近邻高斯过程（见图 A.5，可在线获取）。$\mathscr{S}$ 的另一种选择可以基于树状高斯过程的设置（Gramacy 和 Lee 2008）^[22]。

我们的模拟实验表明: 基于最近邻高斯过程模型的估计和克里金预测，与真正的 Matérn 高斯模型非常相似，即使是缓慢衰减的协方差（参见附录 A8，可在线获取）。 Matérn 协方差函数随距离单调递减，满足 理论筛选条件（theoretical screening conditions？），即能够根据少数邻居进行准确预测（Stein 2002）^[40]。这或许可以解释具有 Matérn 协方差的最近邻高斯过程模型的出色性能。我们还在最近邻的很大一部分与相应位置负相关的场景中，研究了使用波形协方差函数（不满足筛选条件）的最近邻高斯过程模型的性能。最近邻高斯过程估计仍然接近真实模型参数，克里金表面与真实表面非常相似（参见附录 A9，可在线获取）。

大多数波形协方差函数（如阻尼余弦函数或基数正弦函数）产生具有几个小特征值的协方差矩阵。完整的高斯过程模型不能用于此类模型，因为矩阵求逆在数值上不稳定。最近邻高斯过程模型涉及更小的矩阵求逆，并且在某些情况下可实现（例如，对于阻尼余弦模型）。然而，对于基数正弦协方差，最近邻高斯过程也面临数值问题，因为即使是小的 $m×m$ 协方差矩阵，在数值上也是不稳定的。带偏置调整（bias-adjusted）的低秩高斯过程（Finley、Banerjee 和 McRoberts 2009）^[15] 在这方面具有一定的优势，因为协方差矩阵保证具有远离零的特征值，尽管稳定的计算通常需要完整的 Cholesky 分解。

除了可以轻松扩展到具有离散时间的多元和时空设置之外，最近邻高斯过程还可以激发人们在图上基于过程建模的兴趣。示例中包括网络（network)，其中来自结点的数据被假定为与相邻结点相似。在空间数据的区域聚合分析方面，它还提供了新的建模途径和替代马尔可夫随机场模型的方法。此外，当空间和时间被联合建模为使用时空协方差函数的过程时，还有创新的空间。人们需要在空间和时间上构建邻居集，并且探索可扩展性和推断方面的有效策略。

本文方法还需要与其他方法进行比较（例如，参见 Katzfuss 和 Cressie 2012 ^[26]）。最后，正在对替代算法和参数化进行更全面的研究，以实现更快的马尔可夫链蒙特卡洛收敛，包括执行稀疏 Cholesky 分解的直接方法（见第 4 节）。更直接的是，我们计划将较低级别的 C++ 代码迁移到 R 统计环境中的现有 spBayes 包以促进更广泛的用户使用最近邻高斯过程模型。

参考文献

• [1] Banerjee, S., Carlin, BP., Gelfand, AE. Hierarchical Modeling and Analysis for Spatial Data. 2. Boca Raton, FL: Chapman & Hall/CRC; 2014.
• [2] Banerjee S, Finley AO, Waldmann P, Ericcson T. Hierarchical Spatial Process Models for Multiple Traits in Large Genetic Trials. Journal of the American Statistical Association. 2010.
• [3] Banerjee S, Gelfand AE, Finley AO, Sang H. Gaussian Predictive Process Models for Large Spatial Datasets. Journal of the Royal Statistical Society, Series B. 2008.
• [4] Bechtold, WA., Patterson, PL. The Enhanced Forest Inventory and Analysis National Sample Design and Estimation Procedures (SRS-80). Asheville, NC: U.S. Department of Agriculture, Forest Service, Southern Research Station; 2005.
• [5] Bevilacqua M, Gaetan C. Comparing Composite Likelihood Methods Based on Pairs for Spatial Gaussian Random Fields. Statistics and Computing. 2014.
• [6] Crainiceanu CM, Diggle PJ, Rowlingson B. Bivariate Binomial Spatial Modeling of Loa Loa Prevalence in Tropical Africa. Journal of the American Statistical Association. 2008.
• [7] Cressie NAC, Johannesson G. Fixed Rank Kriging for Very Large Data Sets. Journal of the Royal Statistical Society, Series B. 2008.
• [8] Cressie, NAC., Wikle, CK. Statistics for Spatio-Temporal Data. Hoboken, NJ: Wiley; 2011.
• [9] Davis, TA. Direct Methods for Sparse Linear Systems. Philadelphia, PA: Society for Industrial and Applied Mathematics; 2006.
• [10] Diggle PJ, Tawn JA, Moyeed RA. “Model-Based Geostatistics” (with discussion). Applied Statistics. 1998.
• [11] Du J, Zhang H, Mandrekar VS. Fixed-Domain Asymptotic Properties of Tapered Maximum Likelihood Estimators. Annals of Statistics. 2009.
• [12] Eidsvik J, Shaby BA, Reich BJ, Wheeler M, Niemi J. Estimation and Prediction in Spatial Models With Block Composite Likelihoods. Journal of Computational and Graphical Statistics. 2014.
• [13] Emory X. The Kriging Update Equations and Their Application to the Selection of Neighboring Data. Computational Geosciences. 2009.
• [14] Finley AO, Banerjee S, Gelfand AE. spBayes for Large Univariate and Multivariate Point-Referenced Spatio-Temporal Data Models. Journal of Statistical Software. 2015.
• [15] Finley AO, Banerjee S, McRoberts RE. Hierarchical Spatial Models for Predicting Tree Species Assemblages Across Large Domains. Annals of Applied Statistics. 2009.
• [16] Furrer R, Genton MG, Nychka DW. Covariance Tapering for Interpolation of Large Spatial Datasets. Journal of Computational and Graphical Statistics. 2006.
• [17] Gelfand, AE., Banerjee, S. Multivariate Spatial Process Models. In: Gelfand, AE.Diggle, PJ.Fuentes, M., Guttorp, P., editors. Handbook of Spatial Statistics. Boca Raton, FL: Chapman & Hall/CRC; 2010.
• [18] Gelfand AE, Banerjee S, Gamerman D. Spatial Process Modelling for Univariate and Multivariate Dynamic Spatial Data. Environmetrics. 2005.
• [19] Gelfand AE, Ghosh SK. Model Choice: A Minimum Posterior Predictive Loss Approach. Biometrika. 1998.
• [20] Gelfand AE, Kim H-J, Sirmans C, Banerjee S. Spatial Modeling With Spatially Varying Coefficient Processes. Journal of the American Statistical Association. 2003.
• [21] Gramacy, RB., Apley, DW. Local Gaussian Process Approximation for Large Computer Experiments. Available at http://arxiv.org/abs/1303.0383. 2014.
• [22] Gramacy RB, Lee H. Bayesian Treed Gaussian Process Models With an Application to Computer Experiments. Journal of the American Statistical Association. 2008.
• [23] Gramacy, RB., Niemi, J., Weiss, RM. Massively Parallel Approximate Gaussian Process Regression.Available at http://arxiv.org/abs/1310.5182. 2014.
• [24] Higdon, D. Technical Report. Institute of Statistics and Decision Sciences, Duke University; Durham. NC: 2001. Space and Space Time Modeling Using Process Convolutions. 2001.
• [25] Kammann EE, Wand MP. Geoadditive Models. Applied Statistics.52:1–18. 2003.
• [26] Katzfuss M, Cressie N. Bayesian Hierarchical Spatio-Temporal Smoothing for Very Large Datasets. Environmetrics. 2012.
• [27] Kaufman CG, Scheverish MJ, Nychka DW. Covariance Tapering for Likelihood-Based Estimation in Large Spatial Data Sets. Journal of the American Statistical Association. 2008.
• [28] Lauritzen, SL. Graphical Models. Oxford, UK: Clarendon Press; 1996.
• [29] Lin X, Wahba G, Xiang D, Gao F, Klein R, Klein B. Smoothing Spline ANOVA Models for Large Data Sets With Bernoulli Observations and the Randomized GACV. Annals of Statistics. 2000.
• [30] Lindgren, F., Rue, H. and Lindström, J. (2011) ‘An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach’, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4), pp. 423–498. Available at: https://doi.org/10.1111/j.1467-9868.2011.00777.x.
• [31] Moller, J., Waagepetersen, RP. Statistical Inference and Simulation for Spatial Point Processes. 1. Boca Raton, FL: Chapman & Hall/CRC; 2003.
• [32] Rasmussen, CE., Williams, CKI. Gaussian Processes for Machine Learning. 1. Cambridge, MA: The MIT Press; 2005.
• [33] Rue, H., Held, L. Gaussian Markov Random Fields: Theory and Applications. Boca Raton, FL: Chapman & Hall/CRC; 2005.
• [34] Sang H, Huang JZ. A Full Scale Approximation of Covariance Functions for Large Spatial Data Sets. Journal of the Royal Statistical Society, Series B. 2012.
• [35] Schabenberger, O., Gotway, CA. Statistical Methods for Spatial Data Analysis. 1. Boca Raton, FL: Chapman & Hall/CRC; 2004.
• [36] Shaby, BA. The Open-Faced Sandwich Adjustment for MCMC Using Estimating Functions. Available at http://arxiv.org/abs/1204.3687. 2012.
• [37] Shaby BA, Ruppert D. Tapered Covariance: Bayesian Estimation and Asymptotics. Journal of Computational and Graphical Statistics. 2012.
• [38] Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A. Bayesian Measures of Model Complexity and Fit. Journal of the Royal Statistical Society, Series B. 2002.
• [39] Stein, ML. Interpolation of Spatial Data: Some Theory for Kriging. 1. New York: Springer; 1999.
• [40] Stein ML. The Screening Effect in Kriging. Annals of Statistics. 2002.
• [41] Stein ML. Spatial Variation of Total Column Ozone on a Global Scale. Annals of Applied Statistics. 2007.
• [42] Stein ML. A Modeling Approach for Large Spatial Datasets. Journal of the Korean Statistical Society. 2008.
• [43] Stein ML. Limitations on Low Rank Approximations for Covariance Matrices of Spatial Data. Spatial Statistics. 2014.
• [44] Stein ML, Chi Z, Welty LJ. Approximating Likelihoods for Large Spatial Data Sets. Journal of the Royal Statistical Society, Series B. 2004.
• [45] Stroud, JR., Stein, ML., Lysen, S. Bayesian and Maximum Likelihood Estimation for Gaussian Processes on an Incomplete Lattice.Available at http://arxiv.org/abs/1402.4281. 2014.
• [46] Vecchia AV. Estimation and Model Identification for Continuous Spatial Processes. Journal of the Royal Statistical Society, Series B. 1988.
• [47] Vecchia AV. A New Method of Prediction for Spatial Regression Models With Correlated Errors. Journal of the Royal Statistical Society, Series B. 1992.
• [48] Wang Q, Adiku S, Tenhunen J, Granier A. On the Relationship of NDVI with Leaf Area Index in a Deciduous Forest Site. Remote Sensing of Environment. 2005.
• [49] Wikle C, Cressie NAC. A Dimension-Reduced Approach to Space-Time Kalman Fltering. Biometrika. 1999.
• [50] Yeniay O, Goktas A. A Comparison of Partial Least Squares Regression With Other Prediction Methods. Hacettepe Journal of Mathematics and Statistics. 2002.
• [51] Zhang X, Kondraguanta S. Estimating Forest Biomass in the USA Using Generalized Allometric Models and MODIS Land Products. Geophysical Research Letters. 2006.