近似受限似然方法

〖摘 要〗 由于计算负担，似然法通常难以用于大型、位置不规则的空间数据集。即使对于高斯模型，精确计算 $n$ 个观测值的似然也需要 $\mathcal{O}(n^3)$ 运算。任何联合密度都可以写成基于某些观测顺序的条件密度之积，因此一种减少计算的方法是在计算上述条件密度时，仅以部分的 “过去” 观测为条件。本文重点探讨了此类方法如何应用于受限似然的近似，特别展示了如何利用 估计方程方法 判断近似的有效性。此外，过前的工作通常建议以当前观测的历史最近邻观测为条件，但我们通过理论、数值和实例表明，以一些远距离的观测为条件，通常也可以带来相当大的好处。

〖原 文〗 Stein, M.L., Chi, Z. and Welty, L.J. (2004) ‘Approximating likelihoods for large spatial data sets’, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(2), pp. 275–296. Available at: https://doi.org/10.1046/j.1369-7412.2003.05512.x.

1 引言

估计空间协方差结构是地统计学的基础。空间过程经常在不规则的位置点进行观测，这使得对空间过程协方差结构的估计变得相当困难。在地统计学文献中，关于各向同性随机场的二阶结构的推断 通常基于经验变异函数：根据某个点间距离范围内所有观测值对（pairs）计算得到的均方增量，通常表示为距离的函数。由于变异函数必须有条件地负定（Chiles 和 Delfiner，1999）^[2]，实际的变异函数估计通常需要选择变异函数的某种参数化类型，并估计最优的参数值。

（1）高斯场与似然方法

在估计空间协方差函数的参数时，很自然地会考虑基于似然方法（如最大似然）和贝叶斯方法。

对于高斯随机场而言，可以很容易地写出似然函数的精确表达式。尽管基于高斯假设的似然方法尚未被地统计学界完全接受（Chiles 和 Delfiner (1999) ^[2]，第 110 页），但近年来它们已经被更频繁地应用（Lark，2000 年 ^[12]；Pardo-Iguzquiza，1998 年^[16]；Pardo-Iguzquiza 和 Dowd，1997 年^[17]，1998 年^[18]；Park 和 Baek，2001 年^[19]）。此外，对某些非高斯随机场的推断算法来说，计算高斯随机场的似然是必不可少的步骤 (Diggle 等，1998)^[4]。但不幸的是，即使仅限于高斯随机场，当 $n$ 个观测值的位置分布不规则时，似然函数的每次计算都需要 $\mathcal{O}(n^3)$ 次运算，这对于大型数据集来说不切实际。

（2）基于条件分解的近似似然

Vecchia (1988) ^[25] 基于联合密度的乘法法则（即任何联合密度都可以写成条件密度的乘积），建议对空间数据的似然进行简单近似。令符号 $p$ 表示密度（可能是有条件的），假设随机向量 $\mathbf{Z} = (Z_1, \ldots ,Z_n)^{\prime}$ 具有联合密度 $p(\mathbf{Z}；\boldsymbol{\phi})$，其中 $\boldsymbol{\phi}$ 是未知的参数向量。如果将 $\mathbf{Z}$ 划分为（可能不同长度的）子向量 $\mathbf{Z}_1, \ldots ,\mathbf{Z}_b$，并定义向量符号 $\mathbf{Z}_{(j)}= (\mathbf{Z}_1 \ldots \mathbf{Z}_j)$，则根据乘法法则，我们有：

$$p(\mathbf{Z}; \boldsymbol{\phi} ) = p(\mathbf{Z}_1; \boldsymbol{\phi} ) \prod^{b}_{j=2} p(\mathbf{Z}_j \mid \mathbf{Z}_{(j-1)}; \boldsymbol{\phi} ) \tag{1}$$

Vecchia (1988)^[25] 指出，在计算 $p(\mathbf{Z}_j \mid \mathbf{Z}_{(j-1)};\boldsymbol{\phi} )$ 时，没有必要严格地以 $\mathbf{Z}_{(j-1)}$ 的所有分量为条件，而是可以只取部分分量作为条件，这样能够有效减少似然计算中矩阵运算的复杂度。特别地，如果对于 $j = 1, \ldots ,b−1$ 来说，$\mathbf{S}_{(j)}$ 是 $\mathbf{Z}_{(j)}$ 的某个子向量，则联合概率密度可以近似为：

$$p(\mathbf{Z}; \boldsymbol{\phi} ) \approx p(\mathbf{Z}_1; \boldsymbol{\phi} ) \prod^{b}_{j=2} p(\mathbf{Z}_j \mid \mathbf{S}_{(j-1)}; \boldsymbol{\phi} )$$

上式是 Vecchia 近似似然的广义形式， $\mathbf{Z}_j$ 被称为第 $j$ 个预测向量，$\mathbf{S}_{(j-1)}$ 被称为条件向量（或条件集），而 $\mathbf{Z}_{(j)}$ 通常被称为完全条件。

Vecchia (1988) 的原始论文中，只考虑了长度为 $1$ 的预测向量，即 $\mathbf{Z}_j = Z_j$，但其实允许更长的预测向量是有用的。Pardo-Iguzquiza 和 Dowd (1997) ^[17] 描述了一个使用 Vecchia 的方法来近似高斯随机场似然的软件。Eide 等 (2002)^[5] 使用这种形式的近似似然，作为在地震数据贝叶斯分析中实施 MCMC 算法（注：贝叶斯方法也需要计算似然函数）的一个步骤。

注： Vecchia 的原文见 《Vecchia 近似似然法》，但该文献中只考虑了单变量的条件分布（即 $\mathbf{Z}_j = Z_j$），没有考虑子向量的条件分布问题。本文此处对齐进行了泛化。

（3）基础模型设置

在本文整个工作中，我们将假设 $\mathbf{Z} \sim \mathcal{N}\{F\boldsymbol{\beta}, K(\boldsymbol{\theta})\}$，即均值函数为一个线性函数，其中 $F$ 是秩为 $p$ 的已知 $n × p$ 矩阵，$\boldsymbol{\beta} \in \mathbb{R}^p$ 是未知系数的向量；而 $\boldsymbol{\theta} \in \boldsymbol{\Theta}$ 是 $\mathbf{Z}$ 的协方差矩阵的未知参数向量，其长度为 $q$。因此模型的总参数向量为 $\boldsymbol{\phi} = (\boldsymbol{\beta}, \boldsymbol{\theta})$。

注：协方差矩阵的未知参数在二阶平稳假设中，主要指协方差函数的参数，如：

• 常用的各向同性 Matern 协方差函数中，包含边缘方差 $\sigma^2$、变程参数 $\rho$ 或 $\kappa$、平滑参数 $\nu$ 等。
• Vecchia（1988 ^[25]）中采用了用谱密度定义的、支持各向异性的协方差函数，其协方差参数包括：边缘方差 $\sigma^2$、轴缩放因子 $\lambda$、旋转因子 $\alpha$ 等。

如果模型中还包括块金效应（白噪声），则还需要定义噪声方差 $\sigma^2_{\eta}$。

（4）受限最大似然法

为了估计协方差矩阵的参数 $\boldsymbol{\theta}$，最大似然法通常会假设 $\boldsymbol{\beta}$ 已知，但这往往会造成对空间过程变异性的低估（Stein（1999）^[23]， 第 6.4 节 ）。 受限最大似然法 (Restricted maximum likelihood，REML) 在估计 $\boldsymbol{\theta}$ 时，仅使用 ** “比照” （contrast）**，或使用均值不依赖于 $\boldsymbol{\beta}$ 的那些 观测值的线性组合，进而避免了低估可变性的问题。 Kitanidis (1983)^[11] 是第一个建议使用受限最大似然法来估计空间协方差函数参数的人。受限最大似然法和最大似然法的比较将在 第 2 节 的末尾进行。

注： 什么是受限最大似然法？它与最大似然法有何区别？

请参阅 《最大似然法与受限最大似然法的比较》

我们发现 Vecchia (1988)^[25] 的方法也适用于对受限似然方法（高斯观测）的近似。为了实现近似，首先注意，对于高斯场 $\mathbf{Z}$，条件分布 $p(\mathbf{Z}_j \mid \mathbf{Z}_{(j-1)}；\boldsymbol{\phi})$ 作为 $\boldsymbol{\phi}$ 的函数，其作用是依据 $\mathbf{Z}_{(j-1)}$ 来预测 $\mathbf{Z}_j$ 的 最佳线性预测器 (固定效应部分) 的误差的密度。因此，Vecchia 的近似方法其实是用 基于 $\mathbf{S}_{(j-1)}$ 的 $\mathbf{Z}_j$ 最佳线性预测器误差的密度 替换了 基于 $\mathbf{Z}_{(j-1)}$ 的 $\mathbf{Z}_j$ 最佳线性预测器误差的密度。正如可以用最佳线性预测器误差的密度来表示完全似然一样（即概率的乘法法则），受限似然也可以用 最佳线性无偏预测误差的密度 来表示。例如，在给定 $\mathbf{Z}$ 的（ 不包含 $Z_1$ 的 ）子向量 $\mathbf{S}$ 时，$Z_1$ 的最佳线性无偏预测是能够最小化 $Z_1 − \boldsymbol{λ}^{\prime}\mathbf{S}$ 方差的 $\boldsymbol{λ}^{\prime}\mathbf{S}$ 线性组合，其约束条件是：对于 $\boldsymbol{\beta}$ 的所有可能值 $Z_1 − \boldsymbol{λ}^{\prime}\mathbf{S}$ 均值为 $0$。

第 2 节 中的 〖命题 1〗 表明，可以使用类似于 式 (1) 的表达式来计算受限似然。具体来说，如果我们用 $\mathbf{Z}_1$ 的一个线性独立 “比照” 集的联合密度替换 $p(\mathbf{Z}_1； \boldsymbol{\phi} )$，并且用基于 $\mathbf{Z}_{(j-1)}$ 的 $\mathbf{Z}_j$ 最佳线性无偏预测误差的联合密度来替换 $p(\mathbf{Z}_j \mid \mathbf{Z}_{(j-1)}；\boldsymbol{\phi} )$，我们就可以获得受限似然，它仅通过 $\boldsymbol{\theta}$ 依赖于 $\boldsymbol{\phi}$。因此，我们通过 $\mathbf{S}_{(j-1)}$ 而不是 $\mathbf{Z}_{(j-1)}$ 考虑 $\mathbf{Z}_j$ 的最佳线性无偏预测误差密度，进而获得了 Vecchia 近似似然的天然模拟。

（5） 观测数据的排序

为了实现这种近似，我们需要以某种方式对观测结果进行排序，并选择预测向量和条件向量。我们同意 Vecchia (1988) 的观点，即观测数据的排序并不重要，我们使用简单排序（例如观测位置沿某个轴的投影排序）。

Vecchia (1988) 建议 $\mathbf{S}_{(j-1)}$ 由最接近 $Z_j$ 的 $m$ 个观测值组成。这种选择的优点是相当简单。正如 Vecchia (1988) 指出的那样，条件向量 $\mathbf{S}_{(j-1)}$ 的更具统计相关性的选择准则，通常取决于未知参数。这使得很难找到在整个参数空间中都能很好工作的条件向量。不过，在许多情况下，选择包含一些远距离观测值的条件向量，其效率增益可能足够大，以至于值得做出这样的选择。

Jones 和 Zhang (1997)^[9] 讨论了 Vecchia 的方法在时空过程中的应用。由于空间和时间的不可通约性，他们建议在对时空相关函数的一些初步估计基础上，通过相关强度来定义最近邻。这种选择无法解决我们在使用 Vecchia 方法时发现的问题 ( 第 3 节和第 4 节 ），因为我们考虑的相关函数都是距离的单调函数。不过，他们确实指出了当不同坐标不可比时选择条件集合的困难性。我们将在第 6 节的应用中就将面临此问题，其中水平维度的变化与垂直维度的变化存在根本的不同。

（6）采用的评估方法

Lark (2000)^[12] 指出，迄今为止还没有任何尝试来评估 Vecchia 的近似似然对最终参数估计的影响。我们将在第 2 节中展示，估计方程方法提供了一种评估近似似然的天然方法。具体来说，通过将受限近似似然对 $\boldsymbol{\theta}$ 分量的导数设置为 $0$，我们可以得到一组 $\boldsymbol{\theta}$ 的无偏估计方程。

然后，可以使用稳健的信息准则（Heyde (1997)^[8]，第 12 页）来评估各种预测向量和条件向量选择的统计效率。然而，对于足够大的数据集，即使计算此信息度量也会变得令人望而却步，此时，我们建议通过在第 2 节中介绍的抽样方法对其进行近似。

（6）论文结构

第 3 节 考虑了两个简单的例子，其中可以 “比照” 以下三个标准来选择当预测向量是标量时选择给定长度的条件向量：选择与预测量最近的邻居（被预测的数量），选择条件最小化预测误差方差的向量，并选择最适合估计未知参数的条件向量。第一个标准总是独立于参数的真实值给出相同的条件向量，但第二个和第三个通常不会。然而，对于 第 3 节 中的第一个示例，预测误差准则选择相同的条件向量而不考虑真实参数值，并且这种选择不是由最近的邻居组成的。第二个例子表明，预测误差准则可以选择相同的条件向量，而不管参数值如何，但这种条件向量的选择对于参数估计来说是很差的。

第 4 节 通过使用所得估计方程的信息度量来比较选择调节和预测向量的不同方法的效率。所有示例均基于方形区域中 $1000$ 个不规则位置的观测结果。调节向量的长度为 $m = 8、16、32$，并且将其全部、四分之三或一半的分量选择为预测向量的最近邻，而其余分量则为更远的观测值。在许多情况下，选择条件向量的某些分量不是最近邻可以显着提高结果估计器的效率。对于 $m = 8$，有时情况正好相反，但是对于 $m = 32$，选择一些距离较远的观测值几乎一致地优于仅选择最近的邻居。

对于固定协方差函数，当空间相关性最弱时，最近邻设计往往表现更好，但很难进行进一步的概括。在实践中，我们可能希望在最终选择条件向量之前获得初步的参数估计。

第 5 节 讨论了实施我们的程序时的计算问题。我们认为，在某些情况下，两倍于预测集大小的条件集可能是一个不错的选择，这表明具有多个观测值的预测集通常是可取的。

第 6 节 描述了我们的近似受限似然法在密歇根湖超过 $13000$ 次叶绿素水平测量数据集上的应用。数据以不规则的锯齿状模式获取，叶绿素水平在水平和垂直维度上的变化明显不同，这为调节集的选择带来了有趣的挑战。我们发现，在条件集中包含一些远距离观测值可以显着提高某些参数估计的效率。

第 7 节 讨论了这项工作引起的一些计算和推断问题，值得进一步关注。

2 方法论

本节介绍如何使用最佳线性无偏预测将受限似然写成类似于 式 (1) 的形式，然后得出类似于 Vecchia (1988) 的受限近似似然版本。我们表明，关于该近似的未知协方差参数的导数，会产生这些参数的无偏估计方程。此外，估计方程框架提供了一种评估近似效果的天然方法。

2.1 受限似然的定义

我们将在整个工作中假设协方差矩阵 $K(\boldsymbol{\theta})$ 对于所有 $\boldsymbol{\theta} \in \Theta$ 都是正定的。令 $\mathbf{Z}_i$ 的长度为 $n_i$，并取 $F_i$ 为 $F$ 的第 $n_i$ 行，使得 $E(\mathbf{Z}_i) = F_i \boldsymbol{\beta}$。假设 $\text{rank}(F_i) = p$ 并定义 $n_{(j)} = n_1 + \ldots + n_j$。 当 $j>1$ 时，对于所有 $\boldsymbol{\theta} \in \Theta$，随机变量 $\mathbf{Z}_j$ 可以根据 $\mathbf{Z}_{(j−1)}$ 得到最佳线性无偏预测。

当 $j>1$ 时，令 $\boldsymbol{\theta}$ 的函数 $B_j(\boldsymbol{\theta})$ 是一个 $n_j × n$ 的矩阵，使得 $\mathbf{W}_j(\boldsymbol{\theta}) = B_j(\boldsymbol{\theta}) \mathbf{Z}$ 是 $\mathbf{Z}_j$ 基于 $\mathbf{Z}_{(j-1)}$ 的最佳线性无偏预测的误差向量。因此，$B_j(\boldsymbol{\theta})$ 的最后 $n − n_{(j−1)}$ 列等于 $(\mathbf{I}_{n_j} O)$，其中 $\mathbf{I}_{n_j}$ 是 $n_j$ 阶的单位矩阵，$O$ 是 $0$ 矩阵。

当 $j = 1$ 时，将 $B_1(\boldsymbol{\theta})$ 设为大小为 $(n_1 − p) × n$ 且秩为 $n_1 − p$ 的固定矩阵（独立于参数 $\boldsymbol{\theta}$），使得 $\mathbf{W}_1(\boldsymbol{\theta}) = B_1(\boldsymbol{\theta})\mathbf{Z}$ 是 $\mathbf{Z}_1$ 的一组 “比照” 。令 $B(\boldsymbol{\theta})^{\prime} = (B_1(\boldsymbol{\theta})^{\prime} \ldots B_b(\boldsymbol{\theta})^{\prime})$ 并且 $\mathbf{W}(\boldsymbol{\theta})^{\prime} = (\mathbf{W}_1(\boldsymbol{\theta})^{\prime} \ldots \mathbf{W}_b(\boldsymbol{\theta})^{\prime})$。 则 $\mathbf{W}_j(\boldsymbol{\theta}) \sim \mathcal{N}{\mathbf{0}, V_j(\boldsymbol{\theta})}$，其中 $V_j(\boldsymbol{\theta}) = B_j(\boldsymbol{\theta}) K(\boldsymbol{\theta}) B_j(\boldsymbol{\theta})^{\prime}$。

从现在开始，在不引起混淆的情况下，我们将隐藏 $V_j$ 和 $\mathbf{W}_j$ 对 $\boldsymbol{\theta}$ 的依赖。事实证明 $\mathbf{W}_1, \ldots ,\mathbf{W}_b$ 是独立的，导致以下结果（证明见附录 A.1）。

〖命题 1〗 以 $\mathbf{Z}$ 表示的 $\boldsymbol{\theta}$ 的受限对数似然由下式给出

$$\text{rl}(\boldsymbol{\theta}; \mathbf{Z}) = − \frac{n − p}{2} \log(2π) − \frac{1}{2}\sum^{b}_{j=1} \left [\log \{ \det(V_j)\} + \mathbf{W}^{\prime}_j V^{−1}_j \mathbf{W}_j \right]$$

2.2 近似受限似然的定义

现在考虑通过仅根据 $\mathbf{Z}_{(j-1)}$ 的某个子向量计算 $\mathbf{Z}_j$ 的最佳线性无偏预测来近似受限似然，注意确保存在基于该子向量的 $\mathbf{Z}_j$ 的线性无偏预测器。对于 $j>1$，条件向量 $\mathbf{S}_{(j-1)}$ 现在是 $\mathbf{Z}_{(j-1)}$ 的子向量，$\mathbf{Z}_j$ 的最佳线性无偏预测基于该子向量。令 $S$ 为子向量 $\mathbf{S}_{(1)}, \ldots ,\mathbf{S}_{(b−1)}$ 的集合。

定义 $\mathbf{W}_1(S) = \mathbf{W}_1$，对于 $j>1$，$\mathbf{W}_j(S)$ 是 $\mathbf{Z}_j$ 基于 $\mathbf{S}_{(j-1)}$ 的最佳线性无偏预测误差。令 $V_j(S)$ 为 $\mathbf{W}_j(S)$ 的协方差矩阵。考虑 $\text{rl}(\boldsymbol{\theta}；\mathbf{Z})$ 的如下近似形式：

$$\text{rl}(\boldsymbol{\theta}; S) =− \frac{n − p}{2} \log(2π) − \frac{1}{2} \sum^{b}_{j=1} (\log \left [ \det\{ V_j(S) \}\right] + \mathbf{W}_j(S)^{\prime} V_j(S)^{−1} \mathbf{W}_j(S)) \tag{2}$$

与 Vecchia (1988) 相比，上述近似包含两项创新：（1）它适用于受限似然而不是完全似然，后者取决于 $\boldsymbol{\beta}$ 和 $\boldsymbol{\theta}$；（2） Vecchia (1988) 只考虑了每个预测向量只有一个观测值的特殊情况。

2.3 协方差参数的无偏估计

我们可以使用 式（2） 来定义一组 $\boldsymbol{\theta}$ 的无偏估计。记符号 $\partial _k$ 代表 $\partial / \partial \boldsymbol{\theta}_k$，并令 $g_k(S) = \partial _k \text{rl}(\boldsymbol{\theta}; S)$，则我们有：

$$\begin{align*} g_k(S) &= -\frac{1}{2} \sum^{b}_{j=1} \left[ \text{tr}\{V_j(S)^{−1} \partial _kV_j(S)\} + 2 \mathbf{W}_j(S)^{\prime} V_j(S)^{−1} \partial _k \mathbf{W}_j(S) − \mathbf{W}_j(S )^{\prime} V_j(S)^{−1}\{\partial _k V_j(S)\} V_j(S)^{−1} \mathbf{W}_j(S) \right]\\ \mathbf{G}(S) &= (g_1(S) \ldots g_q(S))^{\prime} \tag{3} \end{align*}$$

其中 $\partial _k V_j(S)$ 和 $\partial _k \mathbf{W}_j(S)$ 的明确表达式由『附录 B』 给出。

$\mathbf{G}(S)$ 被称为 $\boldsymbol{\theta}$ 的 估计函数，而 $\mathbf{G}(S) = 0$ 被称为 $\boldsymbol{\theta}$ 的 估计方程。

当 $j>1$ 时，对于 $\Theta$ 中的任意可能的两组参数值 $θ_0$ 和 $θ_1$，残差 $\mathbf{W}_j(θ_0；S) − \mathbf{W}_j(θ_1; S)$ 是仅依赖于子向量 $\mathbf{S}_{(j-1)}$ 的 ** “比照” （contrast）**，因此 $\partial _k \mathbf{W}_j(S)$ 也是 $\mathbf{S}_{(j-1)}$ 的 “比照” ，因此，$\partial _k \mathbf{W}_j(S)$ 和 $\mathbf{W}_j(S)$ 在 $\boldsymbol{\theta}$ 下是独立的，因为最佳线性无偏预测误差独立于观测的所有 “比照” 。此外，$\mathbf{W}_1(S)$ 不依赖于 $\boldsymbol{\theta}$，所以 $\partial _k \mathbf{W}_1(S) = 0$ 并且完全独立于 $\mathbf{W}_1(S)$。从 式 (3) 可以看出，对于所有 $\boldsymbol{\theta} \in \Theta$ 有 $E\{\mathbf{G}(S)\} = 0$，并且 $\mathbf{G}(S) = 0$ 是 $\boldsymbol{\theta}$ 的无偏估计方程。

我们现在可以使用估计方程理论来研究 $\mathbf{G}(S) = 0$ 的解性质，这也为我们提供了一种研究各种 $S$ 子向量有效性的天然方法。如 Heyde (1997) ^[8]所述，将 $\dot{\mathbf{G}}(S)$ 定义为一个 $q × q$ 的矩阵，其第 $(k, l)$ 个元素为 $\partial g_k(S)=\partial θ_l$。估计函数 $\mathbf{G}(S)$ 中关于 $\boldsymbol{\theta}$ 的稳健信息度量为 ^[8]:

$$\mathcal{E}\{\mathbf{G}(S)\} = (E\{\dot{\mathbf{G}}(S)\})^{\prime}\left[ E\{\mathbf{G}(S) \mathbf{G}(S)^{\prime}\} \right]^{−1} E \{\dot{\mathbf{G}}(S) \}$$

估计方程的一个基本目标是使 $\mathcal{E}\{\mathbf{G}(S)\}$ 在非负定矩阵的偏序中变大。如果对 $S$ 没有限制，那么我们可以设置 $\mathbf{S}_{(j)} = \mathbf{Z}_{(j)}$ （ $j = 1, \ldots ,b − 1$ ），在这种情况下 $\mathcal{E}\{\mathbf{G}(S)\}$ 只是 “比照” 的 Fisher 信息矩阵，即偏序下 $\mathcal{E}\{\mathbf{G}(S)\}$ 的最大化器。

当 $n$ 太大时，为了计算 $\text{rl}(θ； \mathbf{Z})$ 时可行，我们可能会寻求在某些约束下将 $\mathcal{E}\{\mathbf{G}(S)\}$ 变大，例如，预测和条件向量的长度受到某些限制。与 Vecchia (1988) 一样，我们寻求产生良好结果的简单规则。 Vecchia 只考虑了长度为 $1$ 且 $m$ 相当小的 $\mathbf{Z}_j$（最多 $m=10$）。当 $j> m$ 时，Vecchia 令 $\mathbf{S}_{(j-1)}$ 是 $\mathbf{Z}_{(j-1)}$ 中 $\mathbf{Z}_j$ 的 $m$ 个最近邻居。 Vecchia (1988) 表明，这种选择在某些情况下能够产生似然的良好近似，但在其示例中，空间相关性仅在远小于观测区域直径的距离处才不可忽略。但 第 4 节 中的示例表明，有时选择条件向量中的一些观测值与预测向量中的观测值相距甚远会有相当大的优势。

现在让我们考虑 $\mathcal{E}\{\mathbf{G}(S)\}$ 的计算。和 $\boldsymbol{\theta}$ 一样，我们隐藏了 $V_j$ 和 $\mathbf{W}_j$ 对 $S$ 的依赖。利用任意 $\mathbf{S}_{(j-1)}$ 的 “比照” 的 $\mathbf{W}_j$ 的独立性，我们再次生成：

$$E\{\partial _l g_k(S) \} = − \sum^{b}_{j=1} \left[\frac{1}{2} \text{tr}\{ V^{-1}_j(\partial _l V_j) V^{-1}_{j} (\partial _k V_j)\} + \text{tr}\{V^{-1}_j \text{cov}(\partial _k\mathbf{W}_j, \partial _l \mathbf{W}^{\prime}_j)\}\right] \tag{4}$$

此外，

$$\begin{align*} E\{g_k(S) g_l(S)\} &= \sum^{b}_{i,j=1} [\text{tr}\{V^{-1}_i \text{cov}(\mathbf{W}_i, \mathbf{W}^{\prime}_j) V^{-1}_j \text{cov}(\partial _l\mathbf{W}_j, \partial _k\mathbf{W}^{\prime}_i)\} \\ &+ \text{tr}\{ V^{-1}_i \text{cov}(\mathbf{W}_i, \partial _l \mathbf{W}^{\prime}_j) V^{-1}_j \text{cov}(\mathbf{W}_j, \partial _k \mathbf{W}^{\prime}_i)\} \\ &− \text{tr} \{V^{-1}_i \text{cov}(\mathbf{W}_i, \mathbf{W}^{\prime}_j) V^{-1}_j (\partial _l V_j) V^{-1}_j \text{cov}(\mathbf{W}_j, \partial _k \mathbf{W}^{\prime}_i)\} \\ &− \text{tr}\{ V^{-1}_j \text{cov}(\mathbf{W}_j, \mathbf{W}^{\prime}_i) V^{-1}_i (\partial _k V_i) V^{-1}_i \text{cov}(\mathbf{W}_i, \partial _l \mathbf{W}^{\prime}_j)\} \\ &+ \frac{1}{2} \text{tr}\{ V^{-1}_i (\partial _k V_i) V^{-1}_i \text{cov}(\mathbf{W}_i, \mathbf{W}^{\prime}_j) V^{-1}_j (\partial _l V_j) V^{-1}_j \text{cov}(\mathbf{W}_j, \mathbf{W}^{\prime}_i)\} ] \end{align*} \tag{5}$$

通过重复应用

$$\begin{align*} \text{cov}(\mathbf{X}^{\prime}_1 A_1 \mathbf{Y}_1, \mathbf{X}^{\prime}_2 A_2 \mathbf{Y}_2) &= \text{tr}\{A_1 \text{cov}(\mathbf{Y}_1, \mathbf{Y}^{\prime}_2) A^{\prime}_2 \text{cov}(\mathbf{X}_2, \mathbf{X}^{\prime}_1)\} \\ &+ \text{tr}\{A_1 \text{cov}(\mathbf{Y}_1, \mathbf{X}^{\prime}_2) A_2 \text{cov}(\mathbf{Y}_2 , \mathbf{X}^{\prime}_1)\} \end{align*}$$

$(\mathbf{X}^{\prime}_1、\mathbf{Y}^{\prime}_1、\mathbf{X}^{\prime}_2、\mathbf{Y}^{\prime}_2)$ 是均值为 $\mathbf{0}$ 的多元高斯，且 $A_1$ 和 $A_2$ 具有适当阶数的固定矩阵。当 $\mathbf{S}_{(j)} = \mathbf{Z}_{(j)}$ （ $j = 1,\dots, b − 1$ ）时，式 (5) 得到大大简化，使得 $\mathbf{G}(S)$ 是得分函数。在这种情况下，当 $i<j$ 时，由于 $\mathbf{W}_i$ 和 $\partial _k\mathbf{W}_i$ 是 “比照” ，并且为 $\mathbf{S}_{(j)}$ 函数，有 $\text{cov}(\mathbf{W}_i, \mathbf{W}^{\prime}_j) = \text{cov}(\partial _k \mathbf{W}_i, \mathbf{W}^{\prime}_j) = O$。因此，$E\{g_k(S) g_l(S)\}$ 简化为 $−E\{\partial _lg_k(S)\}$，并且 $\mathcal{E}\{\mathbf{G}(S)\} = E\{\mathbf{G}(S) \mathbf{G}(S)^{\prime}\}$，即 “比照” 的 Fisher 信息矩阵。

$\mathcal{E}\{\mathbf{G}(S)\}$ 矩阵不仅作为关于估计函数 $\mathbf{G}(S)$ 的信息度量非常有价值，而且对于推断也很有价值。令 $\hat{\boldsymbol{\theta}}(S)$ 是 $\mathbf{G}(S) = 0$ 的解。在许多情况下，对于足够大的样本量，$\hat{\boldsymbol{\theta}}(S)$ 将近似服从均值为 $\boldsymbol{\theta}$ 、协方差矩阵为 $\mathcal{E}\{\mathbf{G}(S)\}^{−1}$ 的高斯（Heyde （1997 年）^[8]，第 4 章）。可用的定理不适用于当前位置不规则的空间数据的情况，但至少可以合理地认为，如果 $\mathcal{E}\{\mathbf{G}(S)\}^{-1}$ 的对角线元素都足够小，对于高斯随机场的观测， $\hat{\boldsymbol{\theta}}(S)$ 应该近似为 $\mathcal{N} (\boldsymbol{\theta}, \mathcal{E}\{\mathbf{G}(S)\}^{-1})$ 。

对于足够大的 $n$，通过 式（3） 精确计算 $E\{\mathbf{G}(S) \mathbf{G}(S)^{\prime}\}$ 将不可行。一种可能的近似是忽略 式 (5) 中的项 $i = j$，因为选择 $S$ 的整个想法是近似设置，其中 对于所有 $j$ 有 $\mathbf{S}_{(j)} = \mathbf{Z}_{(j)}$，在这种情况下 $i = j$ 项全部为 $0$。这相当于 $\mathbf{G}(S)$ 是真实的得分函数。然而，正如我们在第 4.3 节中所述，这可能导致对 $\mathcal{E}\{\mathbf{G}(S)\}$ 过于乐观的近似。

更好的方法是从 $i = j$ 项中随机抽样，并使用抽样理论获得 $E\{\mathbf{G}(S) \mathbf{G}(S)^{\prime}\}$ 的无偏估计量。我们发现以下简单的分层抽样方案通常很有效（请参阅第 4.3 节）。对于每个 $i = 1, \ldots ,b$，无放回地随机选择一些小数值的 $j_1(i), \ldots ,j_r(i)$ outof $1, \ldots ,i − 1, i + 1, \ldots , b$。定义 $u_{kl}(i, j)$ 为 式（5） 中的被加数，对于 $k, l = 1, \ldots, q$，$E\{g_k(S) g_l(S)\}$ 的估计为：

$$\sum^{b}_{i=1} u_{kl}(i, i) + \frac{b-1}{r} \sum^{b}_{i=1} \sum^{r}_{t=1} u_{kl}\{i, j_t(i)\} \tag{6}$$

用于最大似然法的各种预测和条件集选择方法以及估计方程的使用方法同样可以应用于受限最大似然法。

2.4 受限最大似然法的优势

与最大似然相比，受限最大似然有几个优点：

• 首先，在时间序列设置中，理论研究（Kang 等，2003 年 ^[10]）和模拟研究（Wilson，1988 年 ^[28]；McGilchrist，1989 年 ^[15]；Tunnicliffe Wilson，1989 年^[24]）均表明了 受限最大似然法的优越统计特性（尽管需要进一步的工作来验证这些结果是否适用于本文考虑的近似似然）。当 $p$ 很小而 $n$ 很大时，最大似然估计和受限最大似然法估计之间不太可能有太大差异。

• 其次，与精确的受限最大似然法一样，本文方法只需要我们指定 “比照” 的方差。因此，在第 4 节的一些示例和第 6 节的应用中，我们仅指定了变异函数 $\frac{1}{2} \text{var}\{\mathbf{Z}(\mathbf{x}) − \mathbf{Z}(\mathbf{y})\}$。只要均值函数包含未知常数项，变异函数就足以计算近似受限似然。 具体来说，只要公式中需要 $\text{cov}\{\mathbf{Z}(\mathbf{x}), \mathbf{Z}(\mathbf{y})\}$，将其替换为 $− \frac{1}{2} \text{var}\{\mathbf{Z}(\mathbf{x}) − \mathbf{Z}(\mathbf{y})\}$ 即可得到所需的结果。相对而言，如果仅给定变异函数，完全似然和 Vecchia 的近似似然都无法定义。

• 最后，计算（精确或近似）受限最大似然估计在实践中可能会更容易一些，因为完全似然的数值最大化通常需要在估计 $\boldsymbol{\theta}$ 和 $\boldsymbol{\beta}$ 之间进行迭代（Mardia 和 Marshall，1984）^[14]；相比之下，受限最大似然只是在 $\boldsymbol{\theta}$ 上最大化；如果需要，可以根据受限最大似然法对 $\boldsymbol{\theta}$ 的估计结果，显式地估计 $\boldsymbol{\beta}$。

3 最近邻、预测和估计

当所有 $\mathbf{Z}_j$ 的长度均为 $1$ 时，上述方法的关键因素是条件向量 $\mathbf{S}_{(1)}, \ldots,\mathbf{S}_{(b−1)}$ 的选择。我们可能会合理地想象，在 $\mathbf{S}_{(j-1)}$ 的长度约束下，选择使 $V_j(θ、S)$ 尽可能小的向量将是一个好主意。然而，正如 Vecchia (1988)，第 301 页所指出的，该最小化向量通常取决于 $\boldsymbol{\theta}$。正是出于这个原因，Vecchia 建议使用欧几里得距离（或其他距离）最接近待预测位置的那些观测作为向量的选择。此外，即使我们可以选择条件向量来最小化预测方差，这种设计也不一定是估计目的的好选择。

在本节中，我们将描述两个简单示例来演示 “最近邻设计”、“有利于预测的设计”、“有利于估计的设计” 三者之间的潜在冲突。

3.1 示例一

我们的第一个例子表明，为了预测的目的，最近邻设计可能是一致的次优设计。假设 $Z$ 是分数布朗运动：实线上的高斯过程具有未知常量均值，对于所有 $x$ 和 $y$ 实数，$\frac{1}{2} \text{var}{Z(x) − Z(y)} = θ_2 |x − y|^{θ_1}$ 对于某些 $θ_2 > 0$ 且 $θ_1 \in (0, 2)$。我们在 $−1$、$a$ 和 $0$ 处有观测值，其中 $−1<a<0$，我们希望预测 $Z(1)$，但我们只能选择两个观测值来预测 $Z(1)$。

〖命题 2〗 对于 $\boldsymbol{θ} = (θ_1, θ_2)$ 的任何可能值，产生具有最小误差方差的最佳线性无偏预测的长度为 $2$ 的向量是 $(Z(−1), Z(0))$（除非 $θ_1 = 1$，在这种情况下 $(Z(−1),Z(0))$ 和 $(Z(a), Z(0))$ 都将误差方差最小化）。

〖证 明〗

当 $θ_1 = 1$ 时，$Z$ 是布朗运动，结果是基于所有三个观测值的 $Z(1)$ 的最佳线性无偏预测只是 $Z(0)$，因此考虑 $θ_1 \neq 1$。让我们首先证明 $(Z(−1), Z(0))$ 总是优于 $(Z(a), Z(0))$。 $Z(1)$ 基于$Z(0)$ 的最佳线性无偏预测就是 $Z(0)$，其均方误差为 $2θ_2$。如果我们在 $−r$ （ $r > 0$ ）处添加一个的观测值，则最佳线性无偏预测均方误差的相对减少为 $1 − \text{corr}\{Z(1) − Z(0), Z(0) − Z(−r)\}^2$，其中 $\text{corr}$ 表示相关性。定义

$$ρ_{θ_1}(r) = \text{corr}\{Z(1) − Z(0), Z(0) − Z(−r)\} = \frac{1}{2} (r^{1/2} + r^{−1/2})^{θ_1} − \frac{1}{2} (r^{θ_1/2} + r^{−θ_1/2})$$

〖引理 1〗对于除 $θ_1 = 1$ 之外的所有 $θ_1 \in (0, 2)$，$ρ_{θ_1}(r)^2$ 在 $r = 1$ 时达到其唯一最大值，在 $(0, 1)$ 上严格增加，在 $(1, \infty)$ 上严格减少。

证明在 附录 A.2 中给出。由于相对于 $r>0$ 最小化最佳线性无偏预测的均方误差等同于最大化 $ρ_{θ_1}(r)^2$，因此观测 $Z$ 在 $−1$ 和 $0$ 产生的 $Z(1)$ 的均方误差小于观测对于所有可能的 $\boldsymbol{\theta}$，$θ_1 \neq 1$，在 $a$ 和 $0$。

最后，让我们证明当 $θ_1 \neq 1$ 时， $(Z(−1),Z(a))$ 也劣于 $(Z(−1), Z(0))$。对于任何一组观测和预测位置，如果我们将所有相乘点间距离乘以某个正常数 $u$，BLUP 的均方误差乘以 $u^{θ_1}$。设置 $u = 1/(1 − a) < 1$ 从而得出 $(Z(−1), Z(a))$ 不如 $(Z\{−(1 + a)/(1 − a)\}, Z(0))$ 用于预测 $Z(1)$。由于 $−1 < −(1 + a)/(1 − a) < 0$，引理 1 意味着 $(Z\{−(1 + a)/(1 − a)\}，Z(0))$ 又次于 $(Z(−1), Z(0))$，命题 2 成立。

因此，对于具有未知均值和协方差参数 $\boldsymbol{\theta}$ 的分数布朗运动 $Z$，如果我们想通过使用 $Z(0)$ 和 $Z(−r)$ 来预测 $Z(1)$ 对于某些 $r>0$，那么 $r$ 的最佳选择是 $1$ 与 $\boldsymbol{\theta}$ 的值无关。但是，考虑使用受限似然估计基于 $Z(1)$、$Z(0)$ 和 $Z(−r)$ 且 $r>0$ 的 $\boldsymbol{\theta}$。那么 $r = 1$，它给出了三个等间隔的观测值，可能不是一个好的选择，因为等间隔的观测网络在估计协方差结构方面往往效果不佳（Stein（1999）， 第 6.6 节 ，Pettitt 和 McBratney（1993），Lark（2002） ) 和 Zhu (2002))。这实际上似乎是这样。假设我们通过 $I^{11}_r(θ_1)$ 来评估观测网络的质量，$I^{11}_r (θ_1)$ 是受限似然逆 Fisher 信息矩阵的第一个对角线元素的值，值越小表示估计量越好。我们关注 $θ_1$，因为它是更有趣的参数，但如果我们考虑 Fisher 信息逆矩阵的行列式，我们会得到相同的结果，它等于 $θ_2^2 I^{11}_r (θ_1)$。使用 Stein (1999) 第 6 章中的 式 (8) 和冗长但直接的计算，可以证明

$$I^{11}_r (θ_1) = 4\{1 − ρ_{θ_1}(r)^2\}^2 / [\log^2(r) + r^{−θ_1} B \{B − A \log(r)\}]$$

其中 $A = (1 + r)^{θ_1} − 1 − r^{θ_1}$ 和 $B = (1 + r)^{θ_1} \log(1 + r) − r^{θ_1} \log(r)$。将 $I^{11}_r (θ_1)$ 视为各种 $θ_1$ 值的 $r$ 函数的数值评估强烈表明，对于所有 $θ_1 \in (0, 2)$，它通过 $r = 1$ 最大化，尽管我们没有证明。容易证明的是，对于任意固定的 $θ_1 \in (0, 2)$，随着 $r \rightarrow 0$ 说随着 $r \rightarrow \infty$，

$$\log^2(r) I^{11}_r (θ_1)=I^{11}_1 (θ_1) \rightarrow \log^2(2)(1 − 2^{θ_1−2})^{−2}$$

所以对于任何给定的 $θ_1$，选择 $r$ 非常小或非常大都比选择 $r = 1$ 好得多。

3.2 示例二

在上例中，我们选择了观测的一个子集，然后根据它们的精确受限似然估计参数，它不是形 式 (2) 的近似值。接下来让我们考虑一个简单的例子，其中当使用 式 (2) 来近似受限似然时，对于预测来说是一致最优的设计 $S$ 对于估计来说是灾难性的。假设 $Z$ 是实线上的平稳高斯过程，均值未知且 $\text{cov}\{Z(x), Z(y)\} = C K_α (|x − y|)$，其中 $C$ 和 $α$ 均未知，$C>0$ 且 $α \in I$ , 某些开区间。进一步假设，对于每个 $α \in I$，$K_α$ 是 $[0, \infty)$ 上的正且严格递减的函数。对于 $j = 1, \ldots ,n$，我们在 $x_j = Deltaj$ 处观测到 $Z$。我们希望通过使用长度为 414 的预测向量和条件向量来近似 $C$ 和 $α$ 的受限似然，$\mathbf{Z}_1$ 除外，其长度为 $2$。从对 $K_α$ 的假设来看，很明显，在 $Z(\Delta ), \ldots ,Z\{\Delta (j − 1)\}$，对于 $Z(\Delta j)$ 的最佳线性无偏预测给出最小均方误差的是 $Z\{\Delta (j − 1)\}$，与 $C$ 和 $α$ 的值无关。然而，很明显，不可能通过使用这些条件和预测集来估计 $C$ 和 $α$。更具体地说，我们从 式（2） 中得到

$$\text{rl}(θ; S) =− \frac{n − 1 }{2} \log(2π) − \frac{n − 1}{2} \log [2 C\{K_α(0) − K_α(\Delta )\}] − \frac{1}{2C\{K_α(0) − K_α(\Delta )}\} \sum^{n}_{j=2} [Z\{\Delta (j + 1)\} − Z(\Delta j)]^2$$

因此，尽管可以估计乘积 $C\{K_α(0) − K_α(\Delta )\}$，但不可能分别估计 $C$ 和 $K_α(0) − K_α(\Delta )$。如果我们不总是使用最近的过去观测作为条件向量，而是有时使用第二最近的过去观测，那么就有可能同时估计 $C$ 和 $α$。

4 数值结果

4.1 设计和模型

我们可以选择无数的模型、观测网络、预测向量和条件向量来研究 第 2 节 中描述的近似受限似然的一般方法。在这里，我们将只考虑一个观测网络：从 10000 个站点中随机选择的 1000 个站点平面 .i, j/, i, j = 1, \ldots ,100 中的点，然后通过在 [−0:25, 0:25]2 中添加一个随机点来扰动每个点（可以找到确切位置在 http://galton.uchicago.edu/\sim stein/approx-lik.html）。选择这个样本量是因为尽管它相当大，但仍然可以计算出精确受限似然的 Fisher 信息。观测位置不规则，但具有固定的最小观测距离 0.5，这避免了数值困难以及在考虑没有测量误差的模型时彼此非常接近的观测可能出现的统计问题。

我们首先考虑 \mathbf{Z}j 只有一个观测值，j>1 且 \mathbf{Z}1 的长度为 p + 1，因此 \mathbf{Z}{(j)} 的长度为 n(j) = p + j。我们根据观测的坐标之和对观测进行排序，即从观测区域的左下角到右上角。我们采用恒定大小 m 的条件向量（超出最初的几个）。更具体地说，对于 j>m，S(j−p−1)（\mathbf{Z}_j−p = \mathbf{Z}_j 的条件集）由 {x1, \ldots, xj−1} 中 xj 的 m^{\prime} < m 最近邻加上 m − 来自 {x1, \ldots, xj−1} 的 m^{\prime} 个附加点，其到 xj 的排序距离等于 m + l(j − m − 1)=.m − m^{\prime}/ 对于 l = 1, \ldots, m − m ’。取 l = m − m^{\prime}，我们看到我们总是选择最远的“过去”观测。最后，对于 j m，S(j−p−1) 是整个过去。我们将此设计表示为 D(m, m^{\prime})。设计 D(m,0) 在 j>m 的条件向量中仅使用 xj 的 m 个最近邻居。在不考虑方向的情况下按距离排序选择不是最近邻的点会导致条件向量中观测位置的模式相当随意。我们不声称这种方法在任何意义上都接近最优。然而，当我们尝试更系统地选择遥远的过去点的设计时，我们得到了较差的结果。

我们还考虑了一种使用预测向量的方案，该预测向量通过将观测区域划分为包含几乎相等数量的观测值的近似正方形单元而获得。对于具有 J2 n 的正整数 J，首先根据观测值的第一个坐标将 n 个观测值分成 J 个条带，使得每个条带具有近 n=J 个观测值。具体来说，在根据第一个坐标的值分别对观测值进行排序后，对于 j = 1, \ldots, J，条带 j 包含来自 n(j − 1)=J + 0:5 +1 到 nj=J + 0 的观测值：5 ，其中 x 是不大于 x 的最大整数。接下来，将 J 个条带中的每一个分成 J 个矩形，以便每个矩形具有近 n=J2 个观测值。

具体来说，在根据第二个坐标的值对每个条带内的观测值进行重新排序后（这样来自 n(j − 1)=J + 0:5 +1 到 nj=J + 0:5 的观测值仍然构成条带 j），第 j 个条带中的第 i 个矩形包含从 n{.j − 1/J + i − 1}=J2 + 0:5 +1 到 n{.j − 1/J + i}=J2 + 0:5 的观测值。然后条带 j 的第 i 个矩形中的观测值构成预测向量 \mathbf{Z}(J−1)j+i。对于给定 \mathbf{Z}j，将 \mathbf{Z}j 与过去观测值之间的距离定义为该观测值与 \mathbf{Z}j 中观测值之间距离的最小值。对于正整数 m>m^{\prime}，现在基本上像第一个方案一样选择 \mathbf{S}{(j-1)}：如果 n(j−1) 小于 m，则 \mathbf{S}{(j-1)} = \mathbf{Z}{(j-1)}；否则，\mathbf{S}_{(j-1)} 由 m^{\prime} 最近邻点和那些过去的点组成，这些点与 \mathbf{Z}_j 的排名距离等于 m + l(n.j−1) − m/=.m − m^{\prime}/ 对于 l = 1，\ldots , m - m^{\prime}。用 DJ .m, m^{\prime}/ 表示这个设计。在这里，我们考虑 J = 8 的设计，以便将 1000 个观测值划分为 b = 64 个大小为 15 或 16 的子集。

我们将考虑高斯随机场的协方差结构的两个模型。第一个是指数模型，\text{cov}{\mathbf{Z}(\mathbf{x}), \mathbf{Z}(\mathbf{y})} = θ_2 exp(−θ_1|x − y|=θ_2)，第二个是幂律变异函数模型，1 2 \text{var}{\mathbf{Z}(\mathbf{x}) − \mathbf{Z}(\mathbf{y})} = θ_2|x − y|θ_1 。我们大多将均值函数视为未知常数，但我们确实考虑了幂律模型在坐标上呈线性的均值函数。我们还考虑使用方差测量误差 θ3 观测到的幂律模型： 1 2 \text{var}{\mathbf{Z}(\mathbf{x}) − \mathbf{Z}(\mathbf{y})} = θ3 1{|x−y|>0} + θ_2|x − y|θ_1 。除幂律模型外，我们仅报告设计 D(m、m^{\prime}) 的结果。

4.2 相对效率

让我们首先考虑指数模型的结果。此处使用的参数化 θ_2 exp(−θ_1|x − y|=θ_2) 具有 θ_1 描述过程局部变化的特性 (1 2 \text{var}{Z(\mathbf{x}) − Z(\mathbf{y})} \sim θ_1|x − y | as |x − y|\rightarrow0)，而 θ_2 仅在实质上影响与 θ_2=θ_1 相比不小的尺度变化。表 1 给出了 “比照” 的逆 Fisher 信息矩阵的对角线元素与 E{\mathbf{G}(S)}-1 的对角线元素的比率（以百分比表示），用于各种设计和 θ 的值。这些比率是近似似然与完全似然的相对效率的度量，并且由于 Fisher 信息矩阵提供了最大信息，因此它们必须都小于 1。我们在所有情况下都使用 θ_2 = 1 并考虑 θ_1 = 0： 02, 0:1, 0:5, 2，θ_1 的值越小，对应的随机场相关性越强。对于 θ_1，各种设计的效率并不强烈依赖于 θ_1 的真实值或 m^{\prime}=m，即由最近邻组成的条件集的分数。对于 m = 8，所有效率至少为 77%，对于 m = 32，它们都至少为 94%。相反，θ_2 估计的效率在很大程度上取决于 θ_1 的真实值和 m^{\prime}=m。特别是，对于 θ_1 = 0:02 或 θ_1 = 0:1，m^{\prime} = m 的设计远不如 $m^{\prime} <m$ 的设计，对于 $θ_1 = 0:5$ 具有竞争力，而对于 θ_1 = 2 则略胜一筹。

当 θ_1 较小时，观测中关于 θ_2 的信息并不多，尽管它们的数量很多。因此，当 θ_1 = 0:02 时，逆 Fisher 信息矩阵的第二个对角线元素为 0.552，因此即使是 θ_2 的精确受限最大似然法估计量也将是高度可变的。缺少信息是因为 θ_2 对距离小于 θ_2=θ_1（在本例中等于 50）的随机场的波动几乎没有影响。可用的少量信息包含在较长距离的相关性中，因此仅使用最近邻来近似条件密度效果不佳也就不足为奇了。正如 Stein (1999) 所描述的，如果我们只对随机场的空间插值感兴趣，那么当 θ_1 较小时，θ_2 的值几乎没有影响，因此估计 θ_2 时效率的严重损失可能是可以容忍的。然而，为了估计过程的未知均值，更重要的是，为了评估该均值估计值的方差，θ_2 的值至关重要。

指数模型是协方差函数的三参数 Mat ern 模型的子类（Stein，1999）。我们计算了该模型的 D(m、m^{\prime}) 设计的相对效率，并获得了与指数模型定性相似的结果，对于 m^{\prime}=m = 0:75 或 m，所有三个参数的相对效率都更高在我们考虑的所有情况下，^{\prime}=m = 0:5 和 m = 16 或 m = 32。

表 2 给出了幂律变异函数模型 1 2 \text{var}{Z(\mathbf{x}) − Z(\mathbf{y})} = θ_2|x − y|θ_1 和两个不同均值函数的结果。我们首先将均值视为未知常数。对于较小的 0:6 和 1 和 m = 8 的 θ_1 值，m^{\prime} = m 设计是最好的，而对于 m = 16，尤其是$ m = 32$，$m^{\prime} < m$ 的设计具有竞争力甚至更好，尤其是对于 m^{\prime}=m = 0:75。对于较大的 θ_1 值 1.4 和 1.8，结果截然不同。首先，所有 m = 8 的设计都具有较差的相对效率，尤其是对于估计 θ_1。然而，在所有情况下，m^{\prime} = m 的设计显然是最差的，总体而言，m^{\prime}=m = 0:75 的设计略好于 m^{\prime}=m = 0:5 的设计。对于 θ_1 = 1:8，θ_1 和 θ_2 的相对效率之间存在惊人的相似性，而其他 θ_1 值则不存在这种相似性。我们对这些相似之处没有很好的解释，这些相似之处也出现在其他一些表格中。

在试图解释幂律模型的结果之前，有必要检查一下如果均值函数被视为坐标中的未知线性多项式，它们将如何变化。如表 2 所示，在常量均值情况下发现的许多结果现在都颠倒过来了。特别是，$m^{\prime} = m$ 的设计与 $m^{\prime} < m$ 的设计的相对优点现在对于较大的 θ_1 非常好，但对于较小的 θ_1 非常差，特别是关于 θ_1 估计器的相对效率。相同设计的两个不同均值函数的相对效率变化显示出显着差异，具体取决于 m^{\prime}=m.Form^{\prime} = m，当均值函数为线性时，θ_1 = 0:6 或 θ_1 = 1 的相对效率总是较小，最显着的是当 m^{\prime} = m = 8 且 θ_1 = 0:6 时，θ_1 从 84.1% 增加到 38.4%。

相反，当 θ_1 = 1:8 且 m^{\prime} = m = 8 时，θ_1 和 θ_2 的相对效率在从常数函数变为线性均值函数时均从约 25% 增加到 55%。

为了研究长度大于 1 的预测向量的设计效率，我们考虑了 D8(m，m^{\prime})，其中 m = 16 和 m = 32 用于幂律协方差模型和常量均值。我们不考虑 m = 8，因为我们不建议使用比预测向量短的条件向量。 D8(m, m^{\prime}) 设计的性能是否优于 D(m, m^{\prime}) 设计尚不清楚。条件向量的大小在每种情况下都相同，但是由于 D8(m, m^{\prime}) 设计对每个预测向量中的观测值使用了正确的联合条件分布，因此 D8(m, m^{\prime}) 设计可能具有优势。相比之下，条件向量是专门为 D.m 中的每个观测选择的，m^{\prime}/，而对于 D8.m 的预测向量中的所有 15 或 16 个观测，它们必须相同，m^{\prime}/，这可能有利于 D(m , m^{\prime}) 设计。至少在目前的情况下，D(m, m^{\prime}) 和 D8(m, m^{\prime}) 的表现非常相似，所以我们不单独列出 D8(m, m^{\prime}) 的结果。事实上，对于 θ_1 的四个值和 .m 的六个值，m^{\prime}/（m = 16 或 m = 32，m^{\prime}=m = 1, 0:75, 0:5），D8 的相对效率。 m,m^{\prime}/设计比 D(m,m^{\prime})设计不超过 0:01。对于较大的 θ_1、D8(m、m^{\prime}) 通常略胜一筹（相对效率高达 0:1），尤其是对于 m = 16 和 m^{\prime} = 16 或 m^{\prime} = 8。

为了研究测量误差对设计相对效率的影响，我们考虑幂律模型，其幂 θ_1 等于 1 或 1.4，θ_2 = 1，测量误差方差 θ3 等于 0.2 或 1。当 m = 8，m^{\prime} = 8 的设计通常更优，有时甚至更优。当 m = 32 时，m^{\prime} = 32 的设计始终不如 m^{\prime} = 16 或 m^{\prime} = 24 的设计。

这些数值结果显示了一些重要的总体趋势。首先，m 的值越大，$m^{\prime}<m$ 的设计相对于 $m^{\prime} = m$ 的设计往往表现得越好。实际上，对于 m = 32，在此处提供的所有结果中，D(32、32) 优于 D(32、24) 的唯一情况是 θ_1 = 2 的指数模型。对于该模型，相关性消亡得非常快，因此仅对最近的邻居进行调节效果很好也就不足为奇了。然而，即使在这种情况下，$m^{\prime}<m$ 的设计也仅略低于 m^{\prime} = m 的设计，因此对于 θ_1，D(32, 24) 的相对效率与 D(32, 32) 的相对效率在 0.01 以内和θ_2。总的来说，我们发现具有更强空间相关性的模型有利于设计 D.m，$m^{\prime}/ m^{\prime}<m$，尽管具有常数和线性均值函数的幂律模型的不同结果难以解释。

4.3 近似信息

如果 式 (2) 中对受限似然的近似是准确的，我们可能希望用更简单、更容易计算的 E{ ̇ \mathbf{G}(S)} 代替稳健的信息度量 E{\mathbf{G}(S)}，其分量由 式给出(4).然而，在我们检查过的所有情况下，这会产生过于乐观的信息值，因为 E{ ̇ \mathbf{G}(S)}-1 的对角线元素与 E{\mathbf{G}(S)}-1 的相应对角线元素的比率小于 1。这些比率倾向于粗略地跟踪表 1-3 中显示的模式，因此，每当估计效率非常低时，E{ ̇ \mathbf{G}(S)} 往往会严重过度乐观。过度乐观的程度有时可能非常极端：在设计 D(8、8) 当 θ_1 = 1:8 时，在具有常数均值的幂律模型中，该比率为 0:097。因此，我们不建议使用 E{ ̇ \mathbf{G}(S)} 来近似 E{\mathbf{G}(S)}。

如果 E{\mathbf{G}(S)} 的完整计算不可行，我们建议使用表达 式 (6) 给出的采样方法来近似 E{\mathbf{G}(S) \mathbf{G}(S)^{\prime}}，从而近似 E{\mathbf{G}(S)}。我们在表达 式 (6) 中的设计 D(8、8) 和 r = 3 的几个案例中应用了这种方法，发现结果通常是足够的。例如，考虑 θ_1 = 1:8 和常数均值的幂律模型，其中 E{ ̇ \mathbf{G}(S)}-1 低估了 E{\mathbf{G}(S)}-1 的两个对角线元素约 10 倍。我们应用表达 式 (6) 10 次，其中 r = 3 以近似 E{\mathbf{G}(S) \mathbf{G}(S)^{\prime}}

因此 E{\mathbf{G}(S)}，并且在所有 10 种情况下，E{\mathbf{G}(S)}-1 的对角线的近似值都在真实值的 5% 以内。

使用采样来近似 E{\mathbf{G}(S) \mathbf{G}(S)^{\prime}} 的一个吸引人的特性是在实践中很容易检查近似是否有效。具体来说，将 b 块视为采样单元，我们可以计算 E{g_k(S) g_l(S)} 估计值的经验协方差矩阵，其中 k, l = 1, \ldots ,q 由表达式 (6 ).由于 E{ ̇ \mathbf{G}(S)} 是已知的，E{\mathbf{G}(S)}-1 的估计是 E{\mathbf{G}(S) \mathbf{G}(S)^{\prime}} 估计的已知线性变换。因此，可以很容易地获得 E{\mathbf{G}(S)}-1 元素估计值的标准误差。

4.4 预测与估计

第 3 节 中的示例表明，当目标是估计时，选择条件向量以在真实模型下进行良好预测不一定是个好主意。这些目标之间的差异也可能发生在更大的观测网络中。考虑使用与第 4(2 节中相同的包含 1000 个观测值的网络和设计 D.m、m^{\prime})，以便预测向量的长度为 1。由于 V_j 是给定整个过去的最佳线性无偏预测误差的方差，log{V_j.S /=V_j} 是非负的，接近 0 的值表示其误差方差更接近最佳可能的预测。我们将使用这些对数比率的平均值作为我们对特定设计和 θ 值的预测质量的度量。表 4 给出了模型 1 2 \text{var}{Z(\mathbf{x}) − Z(\mathbf{y})} =|x − y|1:4 和常数均值下的这些平均对数比。增加 m 会降低平均对数比，这是预期的，因为使用更长的条件向量应该会产生更接近通过对所有过去的观测进行条件调整而获得的预测方差。更有趣的是，对于较大的 m^{\prime}=m，这些平均值较小，因此在条件向量中选择所有最近的邻居会产生总体上最好的预测。相比之下，表 2 显示当 θ_1 = 1:4 时，m^{\prime}=m = 1 产生的协方差函数参数估计效率最低。

我们看到，至少在定性上，这些结果支持 第 3 节 中的结果。具体来说，在条件向量中包含提供不同空间尺度上随机场信息的观测对于参数估计来说可能是一个好主意，即使它不是一个好主意用于预测。因此，在开发一种选择条件向量以近似似然的直觉时，考虑什么设计会给出好的预测通常是不合适的。

5 计算方面的考虑

本节考虑计算 式 (2) 中的近似似然所需的计算量。我们认为，条件集和预测集的相对大小的合理选择是预测集的大小约为条件集的一半。

为简单起见，让我们假设所有条件集的大小都相同，所有预测集的大小都相同，d = n = b。在近似实现该算法所需的浮点运算时，我们将假设 c + d 很大，但与 n 相比较小。为了进一步简化考虑，我们假设已知 \mathbf{Z} 的平均值为 0，因此 p = 0。当 p 实际上为正时，只要 p 远小于 c + ，所需的额外计算就可以忽略不计 d.

考虑找到第 j 个块对近似似然的贡献所需的计算。定义

var ( \mathbf{S}_{(j-1)} \mathbf{Z}_j ) =K= ( K11 K12 K21 K22 ) ,

抑制协方差对 θ 和 j 的依赖性。我们需要计算最佳线性无偏预测的误差向量 \mathbf{W}_j = \mathbf{Z}j − K21K−1 11 \mathbf{S}{(j-1)}，其协方差矩阵 V_j = K22 − K21K−1 11 K12，二次形式 \mathbf{W}^{\prime}_jV^{-1}_j \mathbf{W}_j 和 det。 V_j/.要获得这些数量，首先计算 K11 的 Cholesky 分解 GG^{\prime}，这需要 13 个 c3 浮点运算（Golub 和 van Loan (1996)，第 144 页）。然后计算 H = G−1K12，这需要 c2d 浮点运算（Golub 和 van Loan (1996)，第 89 页）。计算 \mathbf{W}_j = \mathbf{Z}j − H^{\prime}.G^{\prime}/−1\mathbf{S}{(j-1)} 对整体工作的贡献可以忽略不计。现在 K21K−1 11 K12 = H^{\prime}H 需要利用所得矩阵的对称性进行 cd2 次浮点运算；计算 V_j 是一项微不足道的工作。接下来，计算 V_j 的 Cholesky 分解需要 13 d3 个浮点运算，由此获得 W^{\prime} jV -1 j \mathbf{W}_j 和 det(V_j) 可以忽略不计。总之，所需的浮点运算数为 13 c3 + c2d + cd2 + 13 d3 = 13 .c + d/3 加上低阶项

图 13 .c + d/3 与 K 的 Cholesky 分解所需的相同。实际上，可以证明 \mathbf{W}^{\prime}_jV^{-1}_j \mathbf{W}_j 和 det(V_j) 可以在 13 .c + 中计算 d/3 浮点运算，通过一个过程获得 \mathbf{Z}j 的每个分量的最佳预测器，根据 \mathbf{S}{(j-1)} 和 \mathbf{Z}_j 的先前分量通过顺序更新 Cholesky 分解获得。

由于有 b = n=d 个块，计算 式 (2) 中的近似似然所需的浮点运算总数为.n=3d/.c + d/3。 4.2 节中的结果表明，c 可能比 d 对近似精度更为关键。因此，对于任何给定的 c，选择 d 以最小化所需的计算量是合理的。作为 d 的函数，.n=3d/.c + d/3 通过 d = 1 2 c 最小化，在这种情况下，所需的浮点运算总数为 .9=4/nc2。这与计算完全似然所需的 13 n3 浮点运算相比。对于 第 4 节 中报告的示例，其中 d \approx 16、c = 32 和 n = 1000，我们得到 d \approx 1 2 c、.9=4/nc2 \approx 2:3 × 106 和 13 n3 \approx 3： 3 × 108，所以近似似然需要不到全似然浮点运算的 1%。对于 n 较大的问题，相对节省可能会更加显着。

为了比较选择多个观测的预测集与 Vecchia (1988) 中大小为 1 的预测集的效率，考虑 d = 16 和 c = 32 与 d = 1 和 c = 32。近似浮点运算计数第一个设计是 2304n，第二个是 11 979n，所以无论 n 是多少，d = 1 的设计都需要大约五倍的计算，至少当对每个块分别进行计算时是这样（参见 第 7 节 ）。

6 应用于密歇根湖叶绿素荧光

6.1 介绍

浮游植物是主要存在于海洋和湖泊中的单细胞藻类，是海洋生态系统和全球碳循环的重要组成部分。水样通常用于测量浮游植物的水平，但此类研究中的样本量不可避免地较小，因此提供的有关浮游植物水平空间变化的信息有限。相比之下，叶绿素荧光与浮游植物的水平大致呈线性关系，可以原位测量，并且可以在大空间尺度上以高频（例如每秒观测一次）记录下来。

本例中使用的荧光剖面是 2000 年 3 月中旬在密歇根湖下游盆地获得的，作为情景事件 - 五大湖实验的一部分。荧光计以连续起伏的方式从湖面拖到湖底，然后沿着从近海到湖南端近岸的 25 公里横断面拖回，提供了 13000 多个原位测量值（图 1）。我们按收集时间对观测结果进行排序，当我们提到“最近的邻居”时，我们根据此顺序定义最近的。

6.2 变异函数模型

探索性分析表明没有明显的协变量包含在均值函数中，并且荧光值的对数比原始荧光值更接近高斯分布。使用 h 表示距海岸的距离（以公里为单位），v 表示深度（以米为单位），因此我们将我们的过程 \mathbf{Z}.h，对数（荧光）测量的 v/ 设为高斯分布，具有未知的常量均值和由 a 索引的变异函数参数 θ。 Welty 和 Stein (2003) 描述了导致变异函数模型的探索性分析

γ.h, v/ = θ0 1{h2+v2>0} + θ_1hθ_2 + θ3[1 − Mν {.θ_2 4h2 + θ_2 5 v2 /1=2 }]

其中 Mν.z/ = 21−ν zν Kν .z/=Γ.ν/ 是 Mat ern 相关函数（因此 Mν.0/ = 1），Kν 是 ν 阶的修正贝塞尔函数（Abramowitz 和 Stegun , 1965).包含术语 θ_1hθ_2 以捕捉水平变化，这与垂直变化有根本的不同。

6.3 调节和预测集设计

调节载体。我们将数据分成大小为 d 的连续块，用于预测向量 \mathbf{Z}j。基于 第 5 节 的结果，我们构建了大小为 2d 的条件向量 \mathbf{S}{(j-1)}，是预测向量长度的两倍。

对于 d = 10 和 d = 20，我们考虑了三种选择条件集的方案。我们让条件集完全由最近邻点、最近邻点和远点的混合以及最近邻点、远点和中程点的混合组成。我们让 DN(d)、DNF(d) 和 DNMF(d) 代表这些各自的方案。例如，对应于 DN(10) 的预测集 \mathbf{Z}_j 的条件集由 \mathbf{Z}_j-1 和 \mathbf{Z}_j-2 组成。对于 DNF(10)，它由 \mathbf{Z}_j−1 以及来自 \mathbf{Z}_j−1, \ldots ,\mathbf{Z}1 的 10 个大致等距的块中的一个点组成。对于 DNMF(10)，它由来自 \mathbf{Z}_j−1 的五个点、来自附近具有相似深度的块的五个点（见图 1）和来自过去的 10 个大致等距块的每个点组成。条件和预测集设计的数据和更详细的描述可在 http://galton.uchicago.edu/\sim stein/approx-lik.html 获得。

6.4 参数估计

我们通过使用共轭梯度算法（按等，1992 年）。虽然不能保证近似似然只有一种模式，但我们已经在几对参数（其余三个参数固定）的网格上绘制了表面，没有发现多种模式的迹象。为了简化我们的初始方法，我们没有最大化 ν，Mat ́ ern 协方差函数中的平滑度参数。在考虑了值 0.5、1.0 和 1.5 之后，我们得出结论，根据受限对数似然和 v = 0 附近经验垂直变异函数的形状的比较，v = 1:0 表现得最好。如表达 式 (6) 所示，分层抽样用于估计 E{G(S) \mathbf{G}(S)^{\prime}}，并获得 E{\mathbf{G}(S)}-1 元素的标准误差。

我们发现预测和调节向量设计确实会影响参数估计。鉴于较长向量的计算时间增加，我们建议从具有较小 d 值的设计中获取参数估计值，然后使用这些估计值作为初始值以使用较大的预测集进行最大化。我们首先使用该策略在 DNF(10) 下获得估计值，得到 ^ θNF(10) = .3:49 × 10−4,3:00 × 10−3,1:22, 1:47 × 10−3, 305, 3:28/，然后使用 ^ θNF(10) 作为起始值找到 DNMF(20) 下的估计 ^ θNMF(20) = .3:50 × 10−4,2:34 × 10−3 ,1:06, 1:48 × 10−3, 306, 3:23/.鉴于这些参数描述了水平方向的行为，并且较大设计中的预测条件集包含明显更多的几乎相同深度的点对，因此两种设计下 θ_1 和 θ_2 的估计值差异并不意外

在条件向量中包含一些较远的点并且具有较长预测条件向量的设计会产生 E{\mathbf{G}(S)}-1 的估计值，对角线元素明显较小（表 5），这表明这些设计的参数估计值的可变性较小。当仅使用最近邻点的方案被修改为也包含更远的点时，对角线元素的平方根出现最大下降。当我们考虑 DNF(d ) 或 DNMF(d ) 而不是 DN(d ) 时，对应于描述过程的远程水平行为的参数 θ_1 和 θ_2 的对角线元素显着减少，即使这些设计具有相等长度的条件向量。包括中程点和远点的方案会为某些参数产生较小的对角线元素，尽管这些收益可能会被其他参数的稍大值所抵消。 DNMF(10) 产生比 DNF(10) 更小的对角线元素，但与 DNF(20) 相比，DNMF(20) 产生更小的 θ_1 值，更大的 θ_2 值和相似的其余参数值。

E{\mathbf{G}(S)}-1 的对角线元素的抽样标准误差对于较小的 r 可能很重要。例如，对于 r = 5 和设计 DNMF(20)，我们使用分层随机样本的基本抽样理论计算了 E{\mathbf{G}(S)}-1 对角线元素的标准误差。对于除θ_1 和θ_2 以外的参数，这些标准误差均小于 E{\mathbf{G}(S)}-1 对角线元素的 5%，而θ_1 和θ_2 分别约为 12%和 14%。 E{\mathbf{G}(S)}-1 的对角线元素从一种设计到另一种设计的差异（当两种设计都使用相同的抽样方案时）往往比单个元素本身的变化小得多，这从几乎相同的值中可以明显看出对于 θ_1 和 θ_2 以外的参数，E{\mathbf{G}(S)}-1 的对角线元素，设计 DNF(20) 和 DNMF(20)。

因此，我们建议根据 E{\mathbf{G}(S)}-1 的值差异，使用较小的 r 值来比较设计（具有相同的非对角线采样 j1(i), \ldots ,jr(i)）。一旦确定了一个好的设计（可能有很多），我们建议增加 r 以获得 E{\mathbf{G}(S)}-1 估计中所需的确定性水平。设计最重要的方面似乎是包括远离过去的点；一旦包含了一些较远的点，就可以通过增加预测条件集的大小来进行改进，并且可能在较小程度上通过调整点选择方案以包括中等范围的点来进行改进。

7 讨论

在查看模型、观测位置以及预测和条件集的良好选择之间的关系方面，我们只是触及了皮毛。当然，我们希望找到更少的随机规则来选择条件集中的远距离位置。另一个关键问题是更好地理解测量误差对良好设计的影响。例如，我们期望间隔紧密的观测值更容易估计测量误差方差，这表明在条件集中选择远距离观测值可能有优势，而当测量时，这种聚集可能效率低下误差方差可以忽略不计。

当观测在规则格上时，应该可以获得足够简单设计的渐近效率的解析表达式。我们建议在可行的情况下使用固定域渐近方法 (Stein, 1999)，以反映强烈依赖的相邻观测的情况，我们希望在条件集中使用远距离观测最有帮助。然而，我们不提倡在观测格子上的平稳过程时使用我们的方法，而是建议使用精确方法，利用所得协方差矩阵的块 Toeplitz 结构或光谱近似，例如 Whittle 似然 （Whittle，1954 年）。此外，对于大于几个观测值的调节和预测集，结果通常会非常混乱，而对于非常小的调节和预测集，结果可能无法提供太多洞察力。

设计是否有效在某种程度上取决于计算的完成方式。例如，此处描述的计算假设不同块的计算是分开进行的。当不同块的条件集之间存在大量重叠时，可以利用这一点来减少计算量，例如，只需计算一次两个或多个条件集中的常见观测值的 Cholesky 分解。如果开发了利用这些重叠的算法，那么条件集之间具有更大重叠的设计可能会受到青睐。近似 (2) 可以很容易地在并行处理器上实现。人们可以将预测和调节集分配给不同的处理器，尽管在平衡算术工作与跨处理器传递信息的需要方面会存在有趣的问题。实现的细节可能会影响各种设计的相对计算效率。

最后，我们提到一个有趣的推论问题，该问题是在贝叶斯分析中使用近似似然引起的。在近似值非常好的情况下，我们可以将近似似然视为实际似然，而不会造成太大伤害。然而，特别是如果似然计算只是马尔可夫链蒙特卡洛算法中单个步骤的一部分，我们可能无法提供高度准确的近似，在这种情况下，所得近似后验作为推断工具的地位会令人怀疑。可能可以使用稳健的信息度量以某种方式调整近似似然，但这种方法的细节和证明其有效性提出了相当大的挑战。

参考文献

• [1] Abramowitz, M. and Stegun, I. (1965) Handbook of Mathematical Functions, 9th edn. New York: Dover Publications.
• [2] Chiles, J. and Delfiner, P. (1999) Geostatistics: Modeling Spatial Uncertainty. New York: Wiley.
• [3] Christensen, R. (1996) Plane Answers to Complex Questions: the Theory of Linear Models, 2nd edn. New York: Springer.
• [4] Diggle, P. J., Tawn, J. A. and Moyeed, R. A. (1998) Model-based geostatistics (with discussion). Appl. Statist., 47, 299–350.
• [5] Eide, A. L., Omre, H. and Ursin, B. (2002) Prediction of reservoir variables based on seismic data and well observations. J. Am. Statist. Ass., 97, 18–28.
• [6] Genton, M. G. and Gorsich, D. J. (2002) Nonparametric variogram and covariogram estimation with FourierBessel matrices. Comput. Statist. Data Anal., 41, 47–57.
• [7] Golub, G. H. and van Loan, C. F. (1996) Matrix Computations, 3rd edn. Baltimore: Johns Hopkins University Press.
• [8] Heyde, C. C. (1997) Quasi-likelihood and Its Application: a General Approach to Optimal Parameter Estimation. New York: Springer.
• [9] Jones, R. H. and Zhang, Y. (1997) Models for continuous stationary space-time processes. In Modelling Longitudinal and Spatially Correlated Data (eds T. G. Gregoire, D. R. Brillinger, P. J. Diggle, E. Russek-Cohen, W. G. Warren and R. D. Wolfinger), pp. 289–298. New York: Springer.
• [10] Kang, W., Wan Shin, D. and Lee, Y. (2003) Biases of the restricted maximum likelihood estimators for ARMA processes with polynomial time trend. J. Statist. Planng Inf., 116, 163–176.
• [11] Kitanidis, P. K. (1983) Statistical estimation of polynomial generalized covariance functions and hydrologic applications. Wat. Resour. Res., 19, 909–921.
• [12] Lark, R. M. (2000) Estimating variograms of soil properties by the method-of-moments and maximum likelihood. Eur. J. Soil Sci., 51, 717–728.
• [13] Lark, R. M. (2002) Optimized spatial sampling of soil for estimation of the variogram by maximum likelihood. Geoderma, 105, 49–80.
• [14] Mardia, K. V. and Marshall, R. J. (1984) Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika, 73, 135–146.
• [15] McGilchrist, C. A. (1989) Bias of ML and REML estimators in regression models with ARMA errors. J. Statist. Computn Simuln, 32, 127–136.
• [16] Pardo-Iguzquiza, E. (1998) Maximum likelihood estimation of spatial covariance parameters. Math. Geol., 30, 95–108.
• [17] Pardo-Iguzquiza, E. and Dowd, P. A. (1997) AMLE3D: a computer program for the inference of spatial covariance parameters by approximate maximum likelihood estimation. Comput. Geosci., 23, 793–805.
• [18] Pardo-Iguzquiza, E. and Dowd, P. A. (1998) Maximum likelihood inference of spatial covariance parameters of soil properties. Soil Sci., 163, 212–219.
• [19] Park, J. S. and Baek, J. S. (2001) Efficient computation of maximum likelihood estimators in a spatial linear model with power exponential covariogram. Comput. Geosci., 27, 1–7.
• [20] Pettitt, A. N. and McBratney, A. B. (1993) Sampling designs for estimating spatial variance components. Appl. Statist., 42, 185–209.
• [21] Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W. T. (1992) Numerical Recipes, 2nd edn. Cambridge: Cambridge University Press.
• [22] Shapiro, A. and Botha, J. D. (1991) Variogram fitting with a general class of conditionally non-negative definite function. Comput. Statist. Data Anal., 11, 87–96.
• [23] Stein, M. L. (1999) Interpolation of Spatial Data: Some Theory for Kriging. New York: Springer.
• [24] Tunnicliffe Wilson, G. (1989) On the use of marginal likelihood in time series model estimation. J. R. Statist. Soc. B, 51, 15–27.
• [25] Vecchia, A. V. (1988) Estimation and model identification for continuous spatial processes. J. R. Statist. Soc. B, 50, 297–312.
• [26] Welty, L. J. and Stein, M. L. (2003) Modeling phytoplankton: covariance and variogram model specification for phytoplankton levels in Lake Michigan. In Geostatistics for Environmental Applications. Boston: Kluwer. To be published. (Available from http://galton.uchicago.edu/∼cises/research/cises-tr1.pdf.)
• [27] Whittle, P. (1954) On stationary processes in the plane. Biometrika, 49, 305–314.
• [28] Wilson, P. D. (1988) Maximum likelihood estimation using differences in an autoregressive-1 process. Communs Statist Theory Meth., 17, 17–26.
• [29] Zhu, Z. (2002) Optimal sampling design and parameter estimation of Gaussian random fields. PhD Dissertation. Department of Statistics, University of Chicago, Chicago.