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$ 个样点的坐标。
1.2 问题域及研究情况
（1）问题域
本研究主要考虑 样点不遵循系统性模式的情况（指空间非规则的数据）。
（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)` 导数的递归求解方法，计算修正贝塞尔函数的现成可用程序也很多，此处均不再赘述。
## 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
$$
模型推断的主要目的，就是根据观测数据集获得上述参数的值，其中涉及大量根据数据计算似然的问题，当 $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.
本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来源 西山晴雪的知识笔记！
相关推荐
2023-01-08
SPDE: 高斯场和高斯马尔可夫随机场之间的明确联系
【摘 要】 连续索引的高斯场 (GF) 是空间统计建模和地统计学中最重要的组成部分，通过协方差函数的定义给出了场性质的直观解释。在计算方面，高斯场受到大 $n$ 问题限制，因为密集矩阵的分解计算成本是维度的三次方（$\mathcal{O}(n^3)$）。尽管当前计算能力处于历史最高水平，但这一事实似乎仍然是许多应用中的瓶颈。与高斯场同样中要的，还有一类离散索引的高斯马尔可夫随机场 (GMRF)，其马尔可夫性质导致精度矩阵的稀疏性，从而使我们可以使用稀疏矩阵的数值算法。对于 $\mathbb{R}^2$ 中的场， GMRF 仅使用了一般算法所需时间的平方根（$\mathcal{O}(\sqrt{n^3})$）。 GMRF 由其完整条件分布分布定义，但在这种参数化形势下，其边缘分布性质并不明确。在本文中，我们展示了：对于 Matérn 类型的某些高斯场，（线性）随机偏微分方程的近似随机弱解，可以为 $\mathbb{R}^d$ 上的任何三角形剖分提供在高斯场和 GMRF 之间的显式链接，进而可以将该高斯场表示为基函数的形式。其好处是：我们既能使用高斯场进行建模，又能够利用...
2023-01-08
固定秩克里金法
【摘 要】 非常大的空间数据集的空间统计具有挑战性。数据集的大小 $n$ 会导致计算最优空间预测变量（例如克里金法）出现问题，因为其计算成本为 $n$ 的三次方。此外，大型数据集通常是在大型空间域上定义，因此感兴趣的空间过程通常在该域上表现出非平稳行为。 通过使用一组固定数量的基函数，可以定义一个灵活的非平稳协方差函数族，这产生了我们称为 “固定秩克里金法” 的空间预测方法。具体来说，固定秩克里金法就是此类非平稳协方差函数支撑下的克里金法。当 $n$ 非常大时，它依赖于计算简化，以获得隐空间过程的空间最佳线性无偏预测器及其均方预测误差。基于最小化加权 Frobenius 范数的方法产生协方差函数参数的最佳估计量，然后将其代入固定秩克里金方程。新方法适用于在整个地球上观测到的非常大的臭氧数据集，其中 $n$ 约为数十万。 【原 文】 Cressie, N. and Johannesson, G. (2008) ‘Fixed rank kriging for very large spatial data sets: Fixed Rank Kriging’, Journal...
2023-01-09
近似受限似然方法
〖摘 要〗 由于计算负担，似然法通常难以用于大型、位置不规则的空间数据集。即使对于高斯模型，精确计算 $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....
2023-01-30
Vecchia 近似似然法的通用框架
【摘要】 高斯过程通常用作函数、时间序列和空间场的模型，但它们对大型数据集在计算上不可行。着眼于高斯过程加上加性噪声项的数据建模典型设置，本文提出了 Vecchia (1988) 方法的泛化作为高斯过程近似的框架。我们展示的通用 Vecchia 方法包含了现有许多流行的高斯过程近似特例，并且允许在统一框架内比较不同方法。通过有向无环图模型，我们确定了推断所需矩阵的稀疏性，从而对计算特性有了新的认识。基于这些结果，我们提出了一种新的稀疏通用 Vecchia 近似，它确保了大型空间数据集的计算可行性，但可以产生比原始 Vecchia 方法近似精度更好的结果。文中提供了几个理论结果并进行了数值比较。 【原文】 Katzfuss, M. and Guinness, J. (2021) ‘A general framework for Vecchia approximations of Gaussian processes’, Statistical Science, 36(1). Available at: https://doi.org/10.1214/19-STS755. 1...
2023-02-01
高斯过程预测的 Vecchia 近似
〖摘 要〗 高斯过程 (GP) 是用于地理空间分析、非参数回归和机器学习的高度灵活的函数估计器，但它们在计算上对大型数据集不可行。 高斯过程的 Vecchia 近似已被用于快速估算参数推断的似然。本文研究了在已观测和未观测位置处进行空间预测时的 Vecchia 近似，包括在大型位置集上获得联合预测分布。我们考虑了用于高斯过程预测的通用 Vecchia 框架，其中包含一些新的和已有的特例。我们从理论和数值上研究了这些方法的准确性和计算特性，并且证明了新方法表现出在空间位置总数上的线性计算复杂性。我们表明，框架内的某些选择会对不确定性量化和计算成本产生强烈影响，从而就哪些方法最适合各种设置提出具体建议。我们还将方法应用于叶绿素荧光卫星数据集，表明新方法比现有方法更快或更准确，并削减了预测结果图中不符合实际的伪影。 〖原 文〗 Katzfuss, M. et al. (2020) ‘Vecchia approximations of Gaussian-process predictions’, Journal of Agricultural, Biological and...
2023-02-07
高斯场和非高斯场的随机偏微分方程方法：10 年回顾
【摘 要】 高斯过程和随机场有着悠久的历史，包含了表示空间和时空相关结构的很多方法，例如：协方差函数、谱表示、再生核希尔伯特空间、基于图的模型等。本文介绍了随机偏微分方程方法（SPDE）如何通过 Hilbert 空间投影，将 Matern 协方差模型与其中几种方法建立起联系，并且每种联系在不同情况下都非常有用。除了主要思想的概述之外，本文还讨论了一些重要的扩展、理论、应用和其他新发展。这些方法包括：马尔可夫模型、非马尔可夫模型、非高斯随机场、非平稳场、任意流形上的时空场等，以及实际计算需要考虑的因素。 【原 文】 Lindgren, F., Bolin, D. and Rue, H. (2022) ‘The SPDE approach for Gaussian and non-Gaussian fields: 10 years and still running’, Spatial Statistics, 50, p. 100599. Available at: https://doi.org/10.1016/j.spasta.2022.100599. 1...