Vecchia 近似似然法

〖摘 要〗 介绍了连续域空间过程的参数估计（指均值函数的参数估计）和模型识别（指残差对应的空间过程模型识别）程序。在本文中，空间过程被假定为具有残差的线性模型，且残差服从二阶平稳高斯随机场，同时假定数据由任意采样位置处空间过程的含噪声观测值组成。本文采用了具有椭圆等值线的二维有理密度函数对空间协方差函数进行建模，文中提出的迭代式估计方法可以减轻非格元数据中常规最大似然估计的许多计算困难。

〖原 文〗 Vecchia, A.V. (1988) ‘Estimation and Model Identification for Continuous Spatial Processes’, Journal of the Royal Statistical Society: Series B (Methodological), 50(2), pp. 297–312. Available at: https://doi.org/10.1111/j.2517-6161.1988.tb01729.x.

1 引言

1.1 背景

令 $\{Z(x, y)\}$ 为一个连续域空间过程，其中 $(x, y) \in \mathbb{R}^2$ 代表空间坐标。该过程由包含空间相关误差的线性回归模型控制：

$$Z(x, y) = f^T(x, y) \beta + \xi(x, y) \tag{1}$$

其中 $f^T(x, y)$ 是 $1 \times r$ 的回归函数预测向量，$\beta$ 是 $r \times 1$ 的未知参数向量，$\{\xi(x, y)\}$ 是一个实值二阶平稳空间随机过程，其均值为零，协方差函数如下式所示：

$$\Gamma(u, v) = \text{cov} \left( \xi(x,y), \xi(x + u, Y + v)\right )$$

本文重点关注对协方差函数的参数化模型 $\Gamma(u,v)$ 的识别，以及式（1） 模型中回归参数 $\beta$ 的估计。

用于参数估计和模型识别的观测数据，由 式（1） 的 $n$ 个观测值组成，这些观测值可能包括加性的测量误差 $\eta$ ，即：

$$z_i = Z(x_i, y_i) + \eta_i , \, 1 \leq i \leq n \tag{2}$$

其中假设 $\eta_i$ 是独立同分布的，服从均值为零，方差为 $\sigma_{\eta}^2$ 的高斯分布；$(x_i,y_i)$ 表示第 $i$ 个样点的坐标。

注： 在贝叶斯（或经验）分层模型中，此处的 $z$ 代表数据模型， $Z(x,y)$ 代表过程模型。

1.2 问题域及研究情况

（1）问题域

本研究主要考虑 样点不遵循系统性模式的情况（指空间非规则的数据）。

• 对于具有系统性模式的情况（如格元数据），可以利用 式 (2) 观测值的协方差矩阵来开发更为高效的方法。

• 格元数据也适合采用频率分析方法，但本文介绍的空间过程是基于空间域的分析。

• Ripley (1981)^[17] 很好地回顾了格元数据的空间过程分析。

（2）可能的途径

已有学者提出了非格元数据的统计分析方法：

• Besag (1975)^[1] 提出了基于最近邻法分析非格元数据的统计方法。
• Cliff 和 Ord (1981)^[3] 中描述的离散空间自回归移动平均方案也可以扩展到非格元数据。

但上述方法并不直接对连续空间过程进行建模。

在被采样空间过程代表物理系统的某种随机模型时，通常需要确定连续过程的局部随机特性。例如，在流体力学中，$\{Z(x, y)\}$ 可能代表速度势场，在这种情况下，$\{- \partial Z/ \partial x\}$ 和 $\{- \partial Z/ \partial y\}$ 分别代表 $x$ 和 $y$ 方向的速度场。随机场 $\{Z(x, y)\}$ 在 空间连续体 上的专业特性知识可用于推断速度场的随机特性，但这种能力在离散化模型上有所欠缺，因为离散模型甚至无法定义有效的偏导数。

（3）地统计学方法

最初由 Matheron（1963 年^[14]，1971 年^[15]）开发了用于采矿业的地统计方法，该方法作为 空间相关性估计 和 对连续空间过程插值 的通用方法，已被广泛接受。Joumel 和 Huijbregts （1978 年）提出的泛克里金法 ^[8] 采用了与 式（1） 和 式（2） 大致相同的模型框架，不过 式（2） 中的误差 $\eta_i$ 被称为 块金效应，用于表征小尺度变异，并且 式（1） 中的 $\{\xi(x, y)\}$ 被假定为仅具有平稳增量（即本征平稳假设），其中与增量有关的协方差由从经验数据开发的 变异函数 来建模。

地统计学中用于模型选择和估计的主要工具是 经验变异函数 ^[8]，它基本上是对真实变异函数的一种矩估计，利用最小二乘法将理论变异函数拟合至经验变异函数（如 Cressie（1985）^[4]）。

除了最小二乘方法外，Mardia 和 Marshall (1984) ^[12]、Marshall 和 Mardia (1985) ^[13]、Kitanidis (1983 ^[9], 1987 ^[10]) 和 Stein (1987) ^[18] 还提出了 最大似然估计 、最小范数二次估计 等替代方法。

（4）本文方法

本研究假设 $\xi(x,y)\}$ 为 二阶平稳过程（注：非本征平稳假设），并且可以使用 Vecchia (1985)^[19] 提出的参数化模型来模拟协方差函数 $\Gamma(u, v)$。该参数化模型在对 $\{\xi(x, y)\}$ 的局部协方差结构建模时非常灵活，并且是基于物理而非经验的（此处主要指克里金法的经验方法）。我们将在 第 2 节 中回顾这些模型。

基于 $\{\xi(x,y)\}$ 的高斯过程假设，本文将使用最大似然法。本文的主要贡献在于：

• 提出了一种用于获得最大似然估计的高效近似统计方法（第 3 节 ），即定义了一系列近似的似然函数 $L_m$。当 $m$ 趋近于 $n$ 时，$L_m$ 趋近于常规的似然函数 $L_n$，这种近似似然函数在 $m$ 比较小时，非常容易计算。
• 开发了一种迭代的估计程序。该程序将 $L_1$ 的估计结果用作 $L_2$ 估计的初值，而后者又用作 $L_3$ 估计的初值，依此类推。该迭代过程中每一步计算的统计量，都可以用于确定参数估计的收敛性。除了对收敛性进行评估外，该统计量也可用于在 $\xi(x,y)\}$ 的多个模型之间进行判别，从而从中识别出最佳模型。在第 4 节中，参数估计和模型识别程序被应用于几个合成数据集和一个实际的地下水数据集。
• 介绍了一些用于分析大型非规则格元数据集的可行统计方法及计算细节。不过，在本研究中，计算方面的考虑是次要的，因此未来有望开发出更有效的新型估算程序。

2 协方差函数的参数化模型

本文对协方差函数的参数化建模，采用源自 Vecchia（1985 ^[19]）的各向异性建模方法，其基本思路是在谱域定义协方差函数的各向异性，而后通过逆傅里叶变换转换为时域的协方差函数。

（1）协方差函数的谱域表示

根据 Vecchia（1985 ^[19]），为 $\{ \xi(x,y) \}$ 的空间协方差结构指定的模型，可以用谱密度函数表示为：

$$S(k_1, k_2) = \int_{- \infty}^{+\infty} \int_{- \infty}^{+\infty} e^{iuk_1 + ivk_2} \Gamma(u,v) du dv \tag{3}$$

其中假设协方差函数 $\Gamma(u,v)$ 绝对可积，因此 式 (3) 在普遍意义上存在。谱密度函数的假设形式是

$$S(k_1, k_2) = S(\kappa) = \frac{\sigma^2 \prod^{q}_{j=1} |\kappa^2 + \theta_j|^{2n_j}} {\prod^{p}_{j=1} |\kappa^2 + \phi_j|^{2m_j}} \tag{4}$$

其中

$$\kappa^2 = [ \lambda^{-1}(k_1 \cos \alpha - k_2 \sin \alpha)]^2 + [ \lambda (k_1 \sin \alpha + k_2 \cos \alpha)]^2\tag{5}$$

式中 $p$ 为正整数，$q$ 为非负整数，$n_j$ 和 $m_j$ 为满足 $\boldsymbol{\Sigma}_{m_j} \leq \boldsymbol{\Sigma}_{n_j} + 1$ 的正整数。 $\lambda$ 为轴缩放因子，$\alpha$ 为方向旋转因子，用于参数化各向异性的形态。

式(4) 为二维空间过程有理谱密度函数的一般形式，可以通过坐标系旋转（ $\alpha$ ）和缩放（$\lambda$） 将其转换为二阶各向同性或方向无关的空间过程。该函数类似于一维过程中最终产生自回归移动平均模型的有理谱密度函数（Priestley（1981）^[16]，第 283 页）。Vecchia (1985) ^[19] 考虑了更一般性的情况，其中部分 $\theta_j$ 和 $\phi_j$ 可能是复数。不过，复值参数会带来计算问题，使得估计变得不可行，因此在本研究中， 式（4） 和 式（5） 的参数空间被认为是：

$$\begin{align*} &\theta_j \in \mathbb{R}_1, \\ &\phi_j \in \mathbb{R}_1^+,\\ &\sigma^2 \in \mathbb{R}_1^+,\\ &\lambda \in R_1^+,\\ &\alpha \in [0,\pi/2] \end{align*}\tag{6}$$

其中 $\mathbb{R}_1$ 是实数，$\mathbb{R}_1^+$ 表示正实数。要使 $S(\kappa)$ 成为有效的谱密度函数，$\phi_j$ 必须为正，并且对 $\lambda$ 和 $\alpha$ 的限制不会导致在描述 式 (5) 的椭圆形状时失去一般性。此外，在 式（5） 中假设单个缩放参数 $\lambda$ 并不会失去一般性。一个旋转轴按 $\lambda$ 缩放，另一个轴按 $\lambda^{-1}$ 缩放，从而使得从 $(k_1, k_2)$ 到新坐标系的变换具有单位雅可比。在比较 $\lambda = 1$ 的各向同性模型和 $\lambda \neq 1$ 的各向异性模型时，此性质非常有用。Brewer 和 Mead (1986) 采用了类似的缩放方法对空间相关性建模。

（2）协方差函数的时域表示

与 式（4） 相对应的协方差函数由下式给出：

$$\Gamma(u, v)=\Gamma(r)=\sigma^2(2 \pi)^{-1}(-1)^{M-1} \sum_{j=1}^p\left[\left(2 m_j-1\right) !\right]^{-1} \partial^{2 m_j-1}\left\{w_j K_0\left(r \sqrt{\phi_j}\right)\right\} / \partial \phi_j^{2 m_j-1} \tag{7}$$

其中：

$$\begin{aligned} r^2 & =[\lambda(u \cos \alpha-v \sin \alpha)]^2+\left[\lambda^{-1}(u \sin \alpha+v \cos \alpha)\right]^2 \\ M & =\sum_{j=1}^p 2 m_j \\ w_j & = \frac{\prod_{l=1}^q\left(\theta_l-\phi_j\right)^{2 n_l}} {\prod_{\substack{l=1,l \neq j}}^p\left(\phi_j-\phi_l\right)^{2 m_l}} \end{aligned}$$

$K_0(\cdot)$ 是零阶的第二类修正贝塞尔函数。要计算过程方差 $\Gamma(0)$，$K_0(r\sqrt{\phi_j})$ 需要被替换为 $- \log(\sqrt{\phi_j})$（Vecchia，1985 ^[19]，命题 3）。 Vecchia (1985) 给出了 式(7) 导数的递归求解方法，计算修正贝塞尔函数的现成可用程序也很多，此处均不再赘述。

注：本文中的协方差函数支持可向异性的表达，因此使问题稍微复杂了一些，而作者在文中并没有明确提及这个特点，对于不了解 Vecchia1985 的读者，在第一次阅读时很容易糊涂。

3 参数估计

本文采用最大似然方法进行参数估计，因此必然涉及似然计算的问题。如果基于 式（7） 的各向异性协方差函数，整个数据集的似然计算效率会很低；为此，作者采用了一种近似计算似然的方法，其基本思路是： 利用概率论中联合分布的条件分解法则，将随机场中的多变量联合似然分解成与数据规模相当的条件分布，然后适当地减小条件集的大小，以期提升计算效率。这进一步带来两个方面需要明确的问题：

• 一是随机变量如何排序，完全的条件分解本身没有排序问题，但人为缩小条件集后自然会带来顺序先后的问题，不同排序有可能带来不同的近似效果；
• 二是条件集如何选择的问题，不同的选择准则可能带来不同的近似效果。

下面，我们首先 在 第 3.1 节 看一下如何进行条件分解和条件集选择；然后考虑如何基于近似似然进行参数估计。

3.1 似然函数的近似

假设 式（1） 中的 $\{\xi(x,y)\}$ 是具有 式(4) 形式谱密度的高斯过程，测量误差 $\eta_i$ 独立同分布且服从 $\mathcal{N}(0,\sigma^2_\eta)$，则需要考虑 式(2) 观测值的似然函数计算问题。令数组形式的观测集合由下式给出：

$$z = F \beta + \xi + \eta \tag{8}$$

其中：

• $F$ 是 $n \times r$ 的矩阵，其第 $i$ 行为 $f^T_i$；
• $z^T = [z_1,\ldots,z_n]$；
• $\xi^T = [\xi_1,\ldots,\xi_n]$；
• $\eta^T = [\eta_1,\ldots,\eta_n]$

此时，式(8) 对应的观测序列之间的顺序关系以及采样位置并不重要；我们会在后面进一步讨论对观测值进行排序的有利方式。

（1）似然的参数化形式

在高斯过程假设下，式(8) 对应的似然函数可以参数化为：

$$L(z)=(2 \pi)^{-n / 2}\left(\sigma^2 \gamma_0\right)^{-n / 2}\left|R+\nu^2 I\right|^{-1 / 2} \exp \left[-\left(2 \sigma^2 \gamma_0\right)^{-1} \varepsilon^{\mathrm{T}}\left(R+\nu^2 I\right)^{-1} \varepsilon\right] \tag{9}$$

其中：

$$\begin{aligned} \gamma_0 & =\operatorname{var}\{\xi(x, y)\} / \sigma^2, \\ R & =\operatorname{corr}(\xi) \\ \nu^2 & =\sigma_\eta^2 / \operatorname{var}\{\xi(x, y)\} \end{aligned}$$

且

$$\varepsilon=z-F \beta$$

注意：

上式中 $R= \operatorname{corr}(\xi)$ 指高斯过程的协方差矩阵，$\varepsilon$ 指含噪声（块金）的观测向量，而 $\gamma_0$ 和 $\nu^2$ 均为和超参数有关的标量。

式（9） 中涉及两处矩阵运算，即 $\left|R+\nu^2 I\right|$ 矩阵行列式计算和 $\left(R+\nu^2 I\right)^{-1}$ 的矩阵求逆运算。所谓大 “n” 问题，正是指这两种计算带来的复杂度。当采用迭代法实现最大似然估计或使用贝叶斯方法做参数推断时，都需要重复多次地计算似然，这进一步增加了整体的计算复杂度。如果能够采用某种方式使这两种运算复杂度降低，则会大大提高计算效率。

目前常见的提效方法包括：

• 降维： 以某种准则找到 “$m$” 个（$m \ll n$) 归纳点，使其结果与 $n$ 个点的结果接近；
• 协方差矩阵稀疏化：努力使 $R$ 稀疏化，进而可以采用一些效率更高的稀疏计算方法；
• 精度矩阵稀疏化：精度矩阵是协方差矩阵的逆矩阵，可以考虑直接使其稀疏化来降低上述矩阵求逆运算的复杂度 。

本文提出的方法（以及其他基于近似似然的方法）则稍有不同，

模型推断的主要目的，就是根据观测数据集获得上述参数的值，其中涉及大量根据数据计算似然的问题，当 $n$ 较大时，似然计算效率较低。

（3）似然的简化

要开发近似似然函数，请注意 式（9） 等价于

$$L(z)=\prod_{i=1}^n p\left(z_i \mid \{ z_j, 1 \leqslant j \leqslant i-1 \}\right) \tag{10}$$

其中 $p(\cdot \mid \cdot)$ 表示条件正态概率密度函数。一般来说，如果 $i$ 很大，随机变量集合 $\{z_j, 1 \leq j \leq i-1\}$ 会在预测 $z$ 时包含大量多余和（或）冗余信息。

这意味着 式（10） 中的第 $i$ 个条件密度可能几乎等同于 $p(z_i \mid z_{im})$，其中 $z_{im}$ 是一个向量，由来自 $\{z_j, 1 \leq j \leq i-1\}$ 的若干个（比如说最多 $m$）观测值组成。根据某些准则来选择 $z_{im}$ 中的观测值通常不可行，因为 $z_{im}$ 取决于 式(4) 模型的形式和特定参数值，执行这种准则会导致参数估计的迭代过程严重不稳定，并大大增加计算时间。

（4）选择条件集

选择 $z_{im}$ 的唯一合乎逻辑的、且独立于模型的方法是： 选择在某种意义上位置最接近 $z_i$ 的那些观测值。因此，将 $m$ 阶的近似似然定义为：

$$L_m(z) = \prod^n_{i=1} p(z_i \mid z_{im}) \tag{11}$$

其中 $z_{im}$ 是一个数组，由 $z_1, \ldots, z_{i-1}$ 中在普通欧氏距离意义上最接近 $z_i$ 的 $\min(i - 1, m)$ 个观测值组成，$d_{ij} = \sqrt{[(x_i - x_j)^2 +(y_i - y_j)^2]}，1 \leq j \leq i-1$。$d_{ij}$ 之间的连接可以以任何一致的方式求解。

当 式（5） 中的各向异性参数 $\lambda$ 较大时，可以用 $z_{im}(\lambda, \alpha)$ 来替换 式（11） 中的 $z_{im}$ ，以便获得近似于完全似然的更好估计，其中 $z_{im}(\lambda, \alpha)$ 由来自 $z_1,\ldots,z_{i-1}$ 的、位于以 $(x_i,y_i)$ 为中心的协方差函数（ 式(7)）的最小等值线内的 $\min(i - 1, m)$ 个观测值组成。不过，除非 $\lambda$ 和 $\alpha$ 是固定且已知的，否则此方法的好处将被其灵活性（ 即允许 $z_{im}$ 依赖于 $\lambda$ 和 $\alpha$ 发生变化）所带来的复杂度抵消。

（5）近似似然的计算

将条件正态密度函数 (Graybill, 1976) ^[7] 的公式代入 式（11），并根据相关性表示结果将得到如下近似似然的计算公式：

$$L_m(z)=(2 \pi)^{-n / 2}\left(\sigma^2 \gamma_0\right)^{-n / 2} \prod_{i=1}^n \omega_{i m}^{-1 / 2} \exp \left[-\left(2 \sigma^2 \gamma_0\right)^{-1} \sum_{i=1}^n \omega_{i m}^{-1} e_{i m}^2\right] \tag{12}$$

其中 $\varepsilon_i$、$\gamma_0$ 和 $\nu$ 与 式（9） 中定义相同。而矩阵 $\omega_{i m}$ 和 $e_{i m}$ 定义如下：

$$\begin{aligned} \omega_{i m} & =1+\nu^2-r_{i m}^{\mathrm{T}}\left(R_{i m}+\nu^2 I\right)^{-1} r_{i m}, \\ e_{i m} & =\varepsilon_i-r_{i m}^{\mathrm{T}}\left(R_{i m}+\nu^2 I\right)^{-1} \varepsilon_{i m}, \\ r_{i m}^{\mathrm{T}} & =\operatorname{corr}\left(\xi_i, \xi_{i m}\right) \end{aligned}$$

上式中数组 $\varepsilon_{im}$ 和 $\xi_{im}$ 的定义类似于 $z_{im}$。$R_{i m}$ 为条件集的协方差矩阵，即：

$$R_{i m}=\operatorname{corr}\left(\xi_{i m}, \xi_{i m}\right)$$

可以很容易看出，式（12） 中最复杂的计算是 $\left(R_{i m}+\nu^2 I\right)^{-1}$ 表示的求逆运算，当 $m$ 值很小时，其计算复杂度会大大降低。

（6）排序方法

可以看出，基于概率条件分解的 式（9） 本身对 $z$ 中观测值的排列具有不变性，但 式（12） 却不是，特别是当 $m$ 较小时。因此，我们显然希望有一个系统的排序程序，以使 式（12） 具有唯一的表达式。

从纯计算角度来看，将较小的 $d_{ij}$ 值与较小的 $i - j$ 值相关联是有利的（ $j < i$ ）。为了确定性，我们假设 $z_i$ 相对于 $x_i$ 或 $y_i$ 做了递增排序。数据位置的绘图可以指示哪种排序会使 $d_{ij}$ 和 $i-j$ 之间存在更大的正相关；否则可以选择任一顺序。本文中分析的所有数据集都是按照 $y$ 坐标的递增顺序排列的。

（7）条件集选择原则

式（12） 表示的近似似然函数是一种可以关于参数实现最大化的恰当形式，我们在下一节会讨论其实现方式。但在此之前，我们先考虑一些似然函数的可能替代近似。

• 方案 1: $m$ 值固定，距离阈值不固定。例如，式（11） 中的条件集 $z_{im}$ 可以被替换为最接近 $z_i$ 的 $m$ 个观测值。该方法与 Besag (1975) ^[1] 的伪似然技术一致，$z$ 中观测值的排序无关紧要。然而，式 (12) 具有随着 $m$ 增加而接近真实似然函数的技术优势。
• 方案 2: 距离阈值固定，$m$ 值不固定。另一种方法是让 $z_{im}$ 包含固定距离 $d_m$ 内的所有观测值 $z_j$, $1 \leq j \leq i - 1$，并让 $z_{im}$ 中的观测值的数量发生变化。对于格元数据，这种方法会产生与本文后面相同的结果。然而，对于不规则间隔的数据，基于固定距离方法的近似似然函数会表现出不稳定的波动，这不利于下文提出的迭代估计方法。

3.2 迭代式最大似然估计

对于中等到大规模的观测数据（即 大 “n” 问题），实现 式（9） 中的完全似然关于未知参数的最大化，变得难以实施。 Mardia 和 Marshall (1984) ^[12] 曾在与 式（1） 相同的模型框架内，开发了精确的最大似然估计程序。但他们提到：该程序对于超过 $150$ 点的数据集变得难以实施。因此，需要一种适用于大型数据集的估计程序。此外，式(7) 的各向异性协方差函数的参数化形式，加剧了 式（9） 中完全协方差矩阵 $R$ 的获取复杂度。

本节提出的估计程序可以缓解上述这些困难。 对于固定的 $m$，我们仅需考虑 式（12） 中的近似似然关于未知参数的最大化即可。其中未知参数包括 式(1) 中的 $\beta$、式(4) 中的方差参数 $\sigma^2$、式(5) 中的缩放 $\lambda$、旋转参数 $\alpha$ 、测量误差参数 $\sigma_{\eta}^2$。

如果假设为各向同性模型且观测误差不存在，则有 $\lambda = 1$ 、$\alpha = 0$、 $\sigma^2_{\eta} = 0$。此时，需要估计的参数只有回归系数 $\beta$ 和方差参数 $\sigma^2$。使用微分方法关于参数 $\sigma^2$ 和 $\beta$ 最小化（对数）近似似然 $-2 \log L_m$，就可以确定参数 $\sigma^2_m$ 和 $\beta_m$ 的近似值。

式（12） 对应的对数似然为：

$$-2 \log L_m^*=n[1+\log (2 \pi)]+n \log \left(\sigma_m^2 \gamma_0\right)+\sum_{i=1}^n \log \omega_{i m} \tag{13}$$

通过最小化 式（13）， 可得参数 $\sigma^2_m$ 和 $\beta_m$ 的近似值：

$$\sigma_m^2=\left(n \gamma_0\right)^{-1} \sum_{i=1}^n \omega_{i m}^{-1} \tilde{e}_{i m}^2 \tag{14}$$

$$\beta_m=\left[\sum_{i=1}^n \omega_{i m}^{-1} g_{i m} g_{i m}^{\mathrm{T}}\right]^{-1}\left[\sum_{i=1}^n \omega_{i m}^{-1} g_{i m} h_{i m}\right] \tag{15}$$

其中

$$g_{i m}=f_i-F_{i m}^{\mathbf{T}}\left(R_{i m}+\nu^2 I\right)^{-1} r_{i m} \tag{16}$$

$$h_{i m}=z_i-r_{i m}^{\mathrm{T}}\left(R_{i m}+\nu^2 I\right)^{-1} z_{i m} \tag{17}$$

在 式(14) 中，$\tilde{e}_{im}$ 表示 式（12） 中的 $e_{im}$，但用 $\beta_m$ 代替了其中的 $\beta$；在 式（16） 中，$F_{im}$ 是一个 $\min(i - 1, m) \times r$ 的矩阵，其第 $k$ 行为 $f^T_{k_{im}}$，其中 $k_{im}$ 是 $z_{im}$ 中第 $k$ 个元素对应的索引。

在迭代过程中，有时需要将 $\beta$ 固定为一个不同于 式(15) 中条件最小值的值 $\tilde{\beta}$，例如，采用普通最小二乘估计：

$$\tilde{\beta}_{\mathrm{OLS}}=\left[\sum_{i=1}^n f_i f_i^{\mathrm{T}}\right]^{-1}\left[\sum_{i=1}^n f_i z_i\right] \tag{18}$$

在此情况下，式（13） 和 式（14） 仍然是基本的估计公式，只是其中 $\tilde{e}_{im}$ 代表了 $\beta$ 等于其指定值 $\tilde{\beta}$ 时的 $e_{im}$。

令 式（13） 中的 $d$ 维未知参数集合由下式给出：

$$\underset{(d \times 1)}{\psi}=\left[\theta^{\mathrm{T}}, \phi^{\mathrm{T}}, \lambda, \alpha, \sigma^2_{\eta}\right]^{\mathbf{T}}$$

其中 $\theta$ 是 $q \times 1$ 的向量，$\phi$ 是 $p \times 1$ 的向量，$\lambda$、$\alpha$ 和 $\sigma^2_{\eta}$ 中的某些参数可能会根据模型的指定从 $\psi$ 中省略。式（13） 中的 $\gamma_0$ 是当 $r = 0$ 、$\sigma^2 = 1$ 时，从 式（7） 获得的 $\psi$ 的某些元素的函数。式（13） 关于 $\psi$ 最小化可以得到 $\psi$ 的估计值 $\hat{\psi}_m$ 被称为 $m$ 阶的 近似最大似然估计；与此同时，可以根据 $\hat{\psi}_m$ 从 式（14） 和 式（15） 获得的估计值 $\hat{\beta}_m$ 和 $\hat{\sigma}^2_m$。

在将 $\beta$ 固定为 $\tilde{\beta}$ 时最小化 式（13） 所获得的参数估计，也被表示为 $\hat{\psi}_m$ 和 $\hat{\sigma}^2_m$。为了最小化 式（13） , 由 R. B. Schnabel 编码的基于拟牛顿法解决无约束非线性问题的 Fortran 优化程序 (Dennis and Schnabel, 1983) ^[6] 被用于本文中的所有应用。当 $\beta$ 包含在估计中时， 式（13） 不会直接最小化，而 $\beta_m$ 总是由 式（15） 给出。 也就是说，在 式（13） 和 式 (14) 中，$\beta_m$ 被固定在例如 $\beta_{m0}$ 处，直至找到关于 $\psi$ 的条件最小值，比如 $\psi_{m1}$。然后，在 $\psi_{m1}$ 处计算 式（15） 以获得 $\beta_m$ 的更新值 $\beta_{m1}$，然后将其用于 式（13） 和 式 (14) 以获得 $\psi$ 的新估计 $\psi_{m2}$。这个过程可以重复几次，直到 $\beta_{mk}$ 和 $\beta_{m,k-1}$ 之间的差变小。基于该过程对大数据集的应用，$k$ 值很少需要大于 $2$ 或 $3$。

迭代估计过程从 $m = 1$ 开始，$\beta$ 固定为普通最小二乘估计 $\tilde{\beta}$（式（18））。然后按前述方法计算 $\hat{\psi}_1$。 $\psi$ 的粗略初始值可以通过在几个点处交互式计算 $- 2 \log L^*_1$ 来确定。得到 $\psi_1$ 后，可以将其作为获取 $\psi_2$ 的初始估计，以此类推，直到估计收敛。有助于确定估计何时收敛的统计量是：

$$\Lambda_m = -2 \log [L^*_m(\hat{\psi}_m)] \tag{19}$$

式右侧是从 式（13） 计算得到的，其中 $\beta$ 固定在 $\tilde{\beta}$ 处。随着 $m$ 的增加，$\Lambda_m$ 接近 $- 2 \log \hat{L}$，其中 $\hat{L}$ 是 式(9) 的最大化精确似然函数。我们从大量模拟和实际数据集的分析中得出的经验是，$m$ 的值通常很小，比如 $m^{\prime}$，因此当 $m > m^{\prime}$ 时，$\Lambda_m$ 的波动可以忽略不计。如以下应用所示，从检查中选择 $m^{\prime}$ 通常是一件简单的事情。

在选择 $m^{\prime}$ 的过程中，$\beta$ 一直保持为 $\tilde{\beta}$；此后，$\beta$ 可以被包含在估计中以获得 $\hat{\psi}_{m^{\prime}}$ 和 $\hat{\beta}_{m^{\prime}}$。

此外，只有在选择 $m^{\prime}$ 之后，我们才将各向异性包括在模型中。这会大大节省计算时间，同时通常对选择过程的影响可以忽略不计。如果 $\hat{\lambda}_{m^{\prime}}$ 与单位矩阵有很大不同时，对于某些 $m > m^{\prime}$，可能需要比较 $- 2 \log [L^{*}_{m^{\prime}}(\hat{\psi_{m^{\prime}}})]$ 和 $- 2 \log [L^{*}_{m}(\hat{\psi_{m^{\prime}}})]$ ，看各向异性模型是否需要更大的 $m^{\prime}$ 值。

4 应用

暂略。

5 结束语

本文的参数估计和模型识别程序对于具有几乎任何样本大小的非格元空间数据集在计算上都是可行的，其前提是数据符合 第 1 节 的假设。

本文中最具限制性的假设是：

• 高斯场 $\{ \xi(x , y)\}$ 中的各向异性可以通过坐标旋转和缩放完全定义，如果要研究其他形式的各向异性，可能需要比 式（4） 和 式(5) 更丰富的模型。
• 式（1） 中的 $Z(x, y)$ 均值选择了线性模型。当不希望指定此类函数时，也可以采用一些去除 $Z(x, y)$ 中平稳趋势的方法，例如 Cressie (1986) 中描述的中值抛光方法，然后对去除趋势后的残差结果应用本文提出的方法（此时可以认为 $\beta = 0$ ）。

尽管这些程序主要是根据直觉开发的，但其在合成和实际数据集上的表现始终支持了几个重要的结论：

• 对于大多数采样方案，几乎所有估计和识别所需的信息都包含在 $1 \leq m \leq 10$ 的近似似然 $L_m$ 中;
• 模型识别方法对各向异性的错误指定具有鲁棒性，允许在假设各向同性的情况下选择特定形式的谱密度。
• 迭代估计统计量对于基于点位置处的一组稀疏观测值来识别连续空间过程的平滑特性是高效的。
• 迭代估计程序对于识别数据中的椭圆各向异性是高效的。

参考文献

• [1] Besag, J. (1975) Statistical analysis of non-lattice data. Statistician, 24,179-195.
• [2] Brewer, A. C. and Mead, R. (1986) Continuous second order models of spatial variation with application to the efficiency of field crop experiments. J. R. Statist. Soc. A, 149, 314-348.
• [3] Cliff, A. D. and Ord, J. K. (1981) Spatial Processes, Models and Applications. London: Pion.
• [4] Cressie, N. A. C. (1985) Fitting variogram models by weighted least squares. J. Int. Ass. Math. Geol., 17,563-586.
• [5] Cressie, N. A. C. (1986) Kriging nonstationary data. J. Amer. Statist. Ass., 81, 625--634.
• [6] Dennis, J. E. and Schnabel, R. B. (1983) Numerical Methodsfor Unconstrained Optimization and Nonlinear Equations. Englewood Cliffs: Prentice-Hall.
• [7] Graybill, F. A. (1976) Theory and Application of the Linear Model, p. 106. North Scituate: Duxbury.
• [8] Joumel, A. G. and Huijbregts, C. J. (1978) Mining Geostatistics. London: Academic Press.
• [9] Kitanidis, P. K. (1983) Statistical estimation of polynomial generalized covariance functions and hydrologic applications. Wat. Resour. Res., 19, 909-921.
• [10] Kitanidis, P. K. (1987) Parametric estimation of covariances of regionalized variables. Wat. Resour. Res., 23, 557-567.
• [11] Lenfest, L. W., Jr (1986) Ground-water levels and use of water for irrigation in the Saratoga Valley, south-central Wyoming, 1980--81. Water-Resources Investigation Report 84-4040. Cheyenne: US Geological Survey.
• [12] Mardia, K. V. and Marshall, R. J. (1984) Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika, 71, 135-146.
• [13] Marshall, R. J. and Mardia, K. V. (1985) Minimum norm quadratic estimation of components of spatial covariance. J. Int. Ass. Math. Geol., 17, 517-525.
• [14] Matheron, G. (1963) Principles of geostatistics. Econom. Geol.,58, 1246-1266.
• [15] Matheron, G. (1971) The theory of regionalized variables and its applications. Cahiers du Centre de Morphologie Mathematique, No.5. Fontainbleau: Centre de Morphologie Mathematique.
• [16] Priestley, M. B. (1981) Spectral Analysis and Time Series. London: Academic Press.
• [17] Ripley, B. D. (1981) Spatial Statistics, ch. 5. New York: Wiley.
• [18] Stein, M. L. (1987) Minimum norm quadratic estimation of spatial variograms. J. Amer. Statist. Ass., 82, 765-772.
• [19] Vecchia, A. V. (1985) A general class of models for stationary two-dimensional random processes. Biometrika, 72, 281-291.