使用过程卷积的空间和时空建模
摘要
一个连续的空间模型可以利用一个核或点扩散函数,通过卷积一个非常简单的、可能是独立的过程来构建。这种构建空间过程的方法比通过空间协变函数的规范提供了许多优势。 特别是,这个过程卷积规范导致了计算上的简化,并且很容易扩展到简单的静止模型之外。 本文使用过程卷积模型来构建灵活且能够容纳大量数据的空间和时空模型。 考虑了来自环境监测的数据。
1. 介绍
使用高斯过程对空间数据进行建模是所有地统计分析的共同点。 该领域的一些著名参考文献包括 Matheron (1963)、Journel 和 Huijbregts (1978)、Ripley (1981)、Cressie (1991)、Wackernagel (1995) 和 Stein (1999)。 一种常见的方法是通过协方差函数 c ( ⋅ ) c(·) c(⋅) 对空间依赖性进行建模,因此任何两点之间的协方差仅取决于它们之间的距离。 距离通常是欧几里得,尽管有时会使用其他度量。
在
R
d
R^d
Rd 上创建高斯过程的另一种可选的构造方法是采用 i.i.d.的
R
d
R^d
Rd 在晶格上的高斯随机变量,并将它们与任意内核进行卷积。 图 1 显示了一个使用高斯核对 i.i.d的 高斯噪声进行卷积的简单示例。 在每个维度中将晶格的密度连续增加 2 倍,同时将变量的方差减少
2
d
2^d
2d 倍,从而在
R
d
R^d
Rd 上产生连续高斯白噪声过程。 这个过程的卷积可以使用
R
d
R^d
Rd 中的一些协变函数等效地定义。 虽然通过 i.i.d. 的卷积定义了一个过程。 高斯晶格变量给出了与通过协变函数定义过程非常相似的结果,卷积构造可以很容易地扩展以允许非标准特征,例如非平稳性、边缘效应、降维、非高斯场和替代空间- 时间模型。 这是采用这种方法的主要动机。
2 通过移动平均构建空间模型
可以通过将连续白噪声过程
x
(
s
)
x(s)
x(s),
s
∈
S
s∈S
s∈S 与平滑核
k
(
s
)
k(s)
k(s) 进行卷积,在一般空间(和时间)区域
S
S
S(例如实平面)上构建高斯过程
z
(
s
)
z(s)
z(s) ,因此
z
(
s
)
=
∫
S
k
(
u
−
s
)
x
(
u
)
d
u
,
s
∈
S
z(s)=\int_Sk(u-s)x(u)du,s∈S
z(s)=∫Sk(u−s)x(u)du,s∈S (1)
得到的
z
(
s
)
z(s)
z(s) 的协方差函数仅取决于位移矢量
d
=
s
−
s
′
d = s - s'
d=s−s′ 并由下式给出
c
(
d
)
=
C
o
v
(
z
(
s
)
,
z
(
s
′
)
)
=
∫
S
k
(
u
−
s
)
k
(
u
−
s
′
)
d
u
=
∫
S
k
(
u
−
d
)
k
(
u
)
d
u
c(d)=Cov(z(s),z(s'))=\int_Sk(u-s)k(u-s')du=\int_Sk(u-d)k(u)du
c(d)=Cov(z(s),z(s′))=∫Sk(u−s)k(u−s′)du=∫Sk(u−d)k(u)du ( 2)
作为一种特殊情况,如果
S
S
S 是
R
m
R^m
Rm 并且
k
(
s
)
k(s)
k(s) 是各向同性的,那么
z
(
s
)
z(s)
z(s) 也是各向同性的,其协方差函数
c
(
d
)
c(d)
c(d) 仅取决于
d
d
d 的大小。 在这种情况下,如果
∫
R
p
k
(
s
)
d
s
<
∞
\int_{R^p} k(s)ds < ∞
∫Rpk(s)ds<∞ 且
∫
R
p
k
2
(
s
)
d
s
<
∞
\int_{R^p} k^2(s)ds < ∞
∫Rpk2(s)ds<∞或
c
(
s
)
c(s)
c(s) 是可积且正定的,则平滑核
k
(
d
)
k(d)
k(d) 和协方差函数
c
(
d
)
c(d)
c(d) 之间存在一对一的关系 。 关系是基于傅里叶变换的卷积定理如下所示,
k
(
s
)
→
F
T
K
(
w
)
→
.
2
K
2
(
w
)
→
I
F
T
c
(
s
)
k(s)→^{FT}K(w)→^{.^2}K^2(w)→^{IFT}c(s)
k(s)→FTK(w)→.2K2(w)→IFTc(s)
k
(
s
)
→
I
F
T
C
1
/
2
(
w
)
→
.
C
(
w
)
→
F
T
c
(
s
)
k(s)→^{IFT}C^{1/2}(w)→^{\sqrt.}C(w)→^{FT}c(s)
k(s)→IFTC1/2(w)→.C(w)→FTc(s)
其中FT 和IFT 表示傅里叶变换及其逆变换,函数
⋅
2
·^2
⋅2 和
.
\sqrt.
.· 是逐点应用的。这给出了协变函数
c
(
s
)
c(s)
c(s) 的谱
C
(
ω
)
C(ω)
C(ω) 与其结果核
k
(
s
)
k(s)
k(s) 之间的关系:
C
(
ω
)
C(ω)
C(ω) 是
k
(
s
)
k(s)
k(s) 的傅立叶变换的平方。请注意,如果过程不是各向同性的,则此关系不再是一对一的。在这种情况下,多个内核可以产生相同的协方差函数。 Thiebaux 和 Pedder(1985 年,第 5 章)以及 Barry 和 VerHoef(1996 年)更详细地探讨了移动平均过程和由其变差函数确定的平稳高斯过程之间的对偶性。图 2 显示了为过程
z
(
s
)
z(s)
z(s) 提供标准高斯、指数和球面协变函数的内核。此外,还显示了由 biwieght 核 (Cleveland, 1979) 诱导的协变函数。与高斯核一样,它会产生一个在原点附近相当平坦的协变函数,并且与球面协变函数一样,在固定距离后依赖性会完全消失。此特定内核将用于本文后面的臭氧建模示例。
2.1 处理卷积模型
在这种基于移动平均的方法下,空间过程 z ( s ) z(s) z(s) 的模型由潜在过程 x ( s ) x(s) x(s) 和平滑内核 k ( s ) k(s) k(s) 的模型规范确定。 一个人对潜在过程和平滑内核的选择可以导致许多有趣的建模方法,其中一些在下面给出。
1.非参数协方差建模
移动平均表示的一个吸引力在于可以对平滑核
k
(
s
)
k(s)
k(s) 建模,而不是必须为正定的协变函数
c
(
s
)
c(s)
c(s)。 例如,Barry and Ver Hoef (1996)、Kern (2000) 和 Ver Hoef、Cressie 和 Barry (2000) 为 k(s) 指定了灵活的模型,以构建高斯过程的非标准协方差。 请注意,还有其他方法可以指定灵活类和正定协方差函数——例如,参见 Gelfand 和 Ecker (1997)
。
2.非正态
x
(
s
)
x(s)
x(s)
基本模型的一个简单扩展是修改管理潜在过程
x
(
s
)
x(s)
x(s) 的规范。 例如,Wolpert 和 Ickstadt (1999) 将
x
(
s
)
x(s)
x(s) 指定为 L´evy 过程,而 Hjort (2000) 使用 Strauss 过程。 因此产生的过程
z
(
s
)
z(s)
z(s) 是不正常的。 一个非常简单的例子是平滑 i.i.d 的晶格。 带有平滑核的 Gamma(2, 1) 随机变量,如下所示。 这基本上产生了一个静止的非高斯场。 这种建模方法可能适用于对通常不遵循正态分布的速率或浓度进行建模。
3.限制
x
(
s
)
x(s)
x(s) 的域
可以选择限制
x
(
s
)
x(s)
x(s) 的域,而不是将
x
(
s
)
x(s)
x(s) 指定为连续的白噪声过程。 一个原因可能是为了更好地考虑正在建模的系统中的依赖关系。 例如,Kern (2000) 在模拟切萨皮克湾的氮浓度时使用了这种修改过的过程来解释形状不规则的边缘。 边缘处潜在过程的这种修改也可能伴随着边缘附近区域中内核的修改。
4.降维
也可以选择将潜在模型
x
(
s
)
x(s)
x(s) 限制在位置
s
1
,
.
.
,
s
m
s1, . . , s_m
s1,..,sm 例如在粗晶格上。 在这种情况下,少量参数
x
(
s
1
)
,
.
.
.
,
x
(
s
m
)
x(s_1), . . . , x(s_m)
x(s1),...,x(sm)有效地控制整个空间过程
z
(
s
)
z(s)
z(s),即使
z
(
s
)
z(s)
z(s) 是连续的并且可能需要在数千个位置。 这是我在本文中关注的主要思想。 除了环境监测问题,这种方法在逆问题中也被证明是卓有成效的,其中空间场的经济参数化可以极大地促进计算。
5.非平稳空间协方差
可以扩展基本表示以允许平滑内核随空间位置而变化。 在这种情况下,空间过程由下式给出
z
(
s
)
=
∫
S
k
s
(
u
)
x
(
u
)
d
u
z(s)=\int_Sk_s(u)x(u)du
z(s)=∫Sks(u)x(u)du
其中
k
s
(
u
)
k_s(u)
ks(u) 表示以站点
s
s
s 为中心的内核,其形状也取决于
s
s
s。 通过允许平滑核
k
s
(
u
)
k_s(u)
ks(u)、
s
∈
S
s ∈ S
s∈S 族随空间位置缓慢变化,可以对依赖结构随 s 变化的空间过程进行建模。 在环境应用中使用这种非平稳过程的例子可以在 Higdon (1998), Higdon 等人中找到。 人。 (2000) 和 Hjort 等人。 (2000 年)。
6.时空模型
通过将一般空间
S
S
S 与时间
T
T
T 扩大,以便现在在空间和时间
S
×
T
S × T
S×T 上定义潜在过程和内核,基本公式可以应用于时空模型。 潜在过程是在空间和时间
x
(
s
,
t
)
x(s, t)
x(s,t) 上定义的,平滑内核
k
(
s
,
t
)
k(s, t)
k(s,t) 也是如此。 示例参见 Higdon (1998)。 本文探讨的一个有趣的替代方案是允许潜在过程
x
(
s
,
t
)
x(s, t)
x(s,t) 随时间演变。 除了纯空间核
k
(
s
)
k(s)
k(s) 外,还通过对潜在过程进行空间平滑来构建时空过程。
z
(
s
,
t
)
=
∫
S
k
(
u
−
s
)
x
(
u
,
t
)
d
u
z(s,t)=\int_Sk(u-s)x(u,t)du
z(s,t)=∫Sk(u−s)x(u,t)du (3)
本文的最后一个例子将这个思想与降维结合起来,使得指定的时空模型可以容纳大量数据。
7.构建依赖的空间过程
最后,过程卷积方法提供了一种构建相关空间过程的方法(例如,参见 Ver Hoef 和 Barry (1998) 和 Ver Hoef、Cressie 和 Barry (2000))。 基本思想是构建过程
z
j
(
s
)
z_j(s)
zj(s),它们在构建过程中共享部分公共潜在过程。 例如,两个过程
z
1
(
s
)
z_1(s)
z1(s) 和
z
2
(
s
)
z_2(s)
z2(s) 可以构造为:
z
1
(
s
)
=
∫
S
0
∪
S
1
k
1
(
u
−
s
)
x
(
u
)
d
u
z_1(s)=\int_{S_0∪S_1}k_1(u-s)x(u)du
z1(s)=∫S0∪S1k1(u−s)x(u)du
z
2
(
s
)
=
∫
S
0
∪
S
2
k
2
(
u
−
s
)
x
(
u
)
d
u
z_2(s)=\int_{S_0∪S_2}k_2(u-s)x(u)du
z2(s)=∫S0∪S2k2(u−s)x(u)du
其中底层潜在过程
x
(
s
)
x(s)
x(s) 位于不相交空间
S
=
S
0
∪
S
1
∪
S
2
S = S_0 ∪S_1 ∪ S_2
S=S0∪S1∪S2 的并集上,并且独立于这些单独的子空间。
z
1
(
s
)
z_1(s)
z1(s) 和
z
2
(
s
)
z_2(s)
z2(s) 的依赖性源于它们对
s
∈
S
0
s ∈ S_0
s∈S0 的
x
(
s
)
x(s)
x(s) 的共同依赖性。 下面的图 4 给出了这种方法的示意图。
3.基本空间模型
这里的目标是提出构建灵活且足够易于处理的空间模型的方法,以便可以使用相当大的数据集进行推理。 因此,需要为潜在过程 x ( s ) x(s) x(s) 提供稀疏支持集。 因此,本文不会考虑将 x ( s ) x(s) x(s) 指定为连续白噪声过程的模型。 有关此类示例,请参见 Higdon 等人。 (1999)。
在这里,我们开发了基本模型的公式。 之后会考虑一些修改。 让
y
1
,
.
.
.
,
y
n
y_1, . . ., y_n
y1,...,yn 是在空间位置 S 中
s
1
,
.
.
.
,
s
n
s_1,...,s_n
s1,...,sn 上记录的数据, 也许最简单的空间模型将数据表示为总平均值
μ
μ
μ、空间过程
z
=
(
z
1
,
.
.
.
,
z
n
)
T
z = (z_1, ..., z_n)^T
z=(z1,...,zn)T 和高斯白噪声
ϵ
=
(
ϵ
1
,
.
.
.
,
ϵ
n
)
T
ϵ = (ϵ_1, ... , ϵ_n)T
ϵ=(ϵ1,...,ϵn)T 方差为
σ
ϵ
2
σ^2_ϵ
σϵ2 ,
y
=
μ
+
z
+
ϵ
y=μ+z+ϵ
y=μ+z+ϵ
其中
z
z
z 的元素是空间过程
z
(
s
)
z(s)
z(s) 对数据位置
s
1
,
.
.
.
,
s
n
s_1,...,s_n
s1,...,sn 的限制
我们将
z
(
s
)
z(s)
z(s) 定义为均值零高斯过程。 但不是通过其协方差函数指定
z
(
s
)
z(s)
z(s),而是由潜在过程
x
(
s
)
x(s)
x(s) 和平滑核
k
(
s
)
k(s)
k(s) 确定。 我们将潜在过程
x
(
s
)
x(s)
x(s) 在空间位置
ω
1
,
.
.
.
,
ω
m
ω_1,..., ω_m
ω1,...,ωm处限制为非零, 也在 S 中并定义
x
=
(
x
1
,
.
.
.
,
x
m
)
T
x = (x_1, . . . , x_m)^T
x=(x1,...,xm)T 其中
x
j
=
x
(
ω
j
)
,
j
=
1
,
.
.
.
,
m
x_j = x(ω_j ), j = 1, . . . ,m
xj=x(ωj),j=1,...,m。 然后将每个
x
j
x_j
xj 建模为来自
N
(
0
,
σ
x
2
)
N(0, σ^2_x)
N(0,σx2) 分布的独立抽取。 然后得到的连续高斯过程是
z
(
s
)
=
∑
j
=
1
m
x
j
k
(
s
−
ω
j
)
z(s) =\sum_{j=1}^m x_jk(s − ω_j)
z(s)=∑j=1mxjk(s−ωj) (4)
其中
k
(
⋅
−
ω
j
)
k(· − ω_j)
k(⋅−ωj) 是以
ω
j
ω_j
ωj 为中心的核。 对于本文中考虑的应用,平滑内核
k
(
⋅
)
k(·)
k(⋅) 将是径向对称内核,例如二元高斯密度或图 2 中的任何其他内核。
这给出了线性模型
y
=
μ
1
n
+
K
x
+
ϵ
y = μ1_n + Kx + ϵ
y=μ1n+Kx+ϵ
其中
1
n
1_n
1n 是 1 的
n
n
n 向量,
K
K
K 的元素由下式给出
K
i
j
=
k
(
s
i
−
ω
j
)
x
j
,
K_{ij} = k(s_i − ω_j)x_j ,
Kij=k(si−ωj)xj,
x
x
x~
N
(
0
,
σ
x
2
I
m
)
,
N(0, σ^2_xI_m),
N(0,σx2Im),
and
ϵ
ϵ
ϵ ~
N
(
0
,
σ
ϵ
2
I
n
)
N(0, σ^2_ϵI_n)
N(0,σϵ2In)
这是一个基本的混合效应模型。 可以使用统计包进行推理,该统计包使用基于可能性的方法用于一般混合模型,例如 SAS 的 proc mixed(Wolfinger 等人,1996)或 Splus 中的 lme(Pinhero 和 Bates,1999)。 请参阅附录中的 Splus 5.0 代码以在此处拟合一维模型。
基于直接似然的方法的替代方法是贝叶斯方法。 完全贝叶斯方法需要对剩余参数 μ μ μ、 σ x 2 σ^2_x σx2 和 σ ϵ 2 σ^2_ϵ σϵ2 进行先验规范。 “默认”公式会在 μ ( π ( μ ) ∝ 1 ) μ (π(μ) ∝1) μ(π(μ)∝1) 之前给出不正确的均匀分布,而在 σ x − 2 σ^{−2}_x σx−2 和 σ ϵ − 2 σ^{−2}_ϵ σϵ−2 之前给出相当平坦的 gamma。 在此公式下没有可用的封闭形式解决方案,因此将需要 MCMC 或其他一些方法来探索所得的后验分布。 如果固定比率 σ x / σ ϵ σ_x/σ_ϵ σx/σϵ,则可以以封闭形式获得 x x x 的后验分布(多元 t 分布)。 对于数据丰富的应用程序,预先估计 σ x / σ ϵ σ_x/σ_ϵ σx/σϵ 然后将其视为固定值可能是相当合理的。
3.1 简单的一维示例
考虑数据对 : ( y i , s i ) (y_i, s_i) (yi,si), i = 1 , . . , n = 30 i = 1, . . , n = 30 i=1,..,n=30 如图 5 所示。通过将 k ( s ) k(s) k(s) 定义为标准差为 2 的单变量高斯分布,并定义潜在过程支持使得 ω j s ω_js ωjs 为 m = 7 m =7 m=7,构建了过程卷积模型 (5) 7 个等距点,范围从 -1 到 12(如图 5 的下部所示。
请注意, ω j s ω_js ωjs 的内核宽度和间距的这种组合通过 (4) 产生了一个空间过程 z ( s ) z(s) z(s),该过程几乎是静止的。 如果间距变得更大,或者如果内核宽度减小,则 z ( s ) z(s) z(s) 的协方差结构会受到稀疏伪影的过度影响
使用 REML (Patterson and Thompson, 1971) 拟合得到的混合效应模型,拟合值(最佳线性无偏预测变量)由图顶部的实线给出。 该图的底部显示了
ω
j
s
ω_js
ωjs(黑点)的位置,还显示了每个
x
j
xj
xj 的估计值。 从(5)拟合值可以表示为
y
′
=
∑
j
=
1
m
K
j
x
j
y'=\sum_{j=1}^mK^jx_j
y′=∑j=1mKjxj
其中
K
j
K^j
Kj 是等式 (5) 的矩阵 K 的第 j 列。 图 5 底部的线条显示了项
K
j
x
j
,
j
=
1
,
.
.
,
7.
K^jx_j , j = 1, . . , 7.
Kjxj,j=1,..,7.
3.2 二维示例
作为一个更严重的例子,我们考虑 1999 年春末一天美国东部的最大 8 小时平均臭氧浓度(图 6)。在查看了数据的一些经验变异函数后,我决定将
k
(
s
)
k(s)
k(s) 指定为范围为 9 度的各向同性二维三角核。因此,诱导的协方差函数在大约 15 度处消失。选择平滑核的更详细方法可以在 Barry 和 VerHoef (1996)、VerHoef 和 Barry (1998) 以及 VerHoef、Cressie 和 Barry (2000) 中找到。选定的
ω
j
s
ω_js
ωjs 由图 6 右侧框架中的黑点显示。它们排列在各向同性六边形网格上,因此每个内部
ω
j
ω_j
ωj 与其他三个
ω
j
′
s
ω_{j'} s
ωj′s 等距。再次通过 REML 获得拟合表面。请注意,当
x
(
s
)
x(s)
x(s) 的支持像这里一样稀疏分布时,直线最小二乘拟合几乎与 REML 拟合相同。这些拟合值的二阶属性尚未进行比较。
4.多分辨率模型
可以通过将空间过程
z
(
s
)
z(s)
z(s) 表示为组件过程的总和来制定多分辨率模型
z
(
s
)
=
∑
ℓ
=
1
p
z
ℓ
(
s
)
z(s)=\sum^p_{ℓ=1}z_ℓ(s)
z(s)=∑ℓ=1pzℓ(s)
其中每个
z
ℓ
(
s
)
z_ℓ(s)
zℓ(s) 都通过其自己的过程卷积模型表示,该模型由对 {
x
ℓ
(
s
)
,
k
ℓ
(
s
)
x_ℓ(s), k_ℓ(s)
xℓ(s),kℓ(s)},
ℓ
=
1
,
.
.
.
,
p
ℓ = 1, ... ,p
ℓ=1,...,p 每个潜在过程
x
ℓ
(
s
)
x_ℓ(s)
xℓ(s) 都支持空间站点
ω
ℓ
1
,
.
.
,
ω
ℓ
m
ℓ
ω_{ℓ1}, . . , ω_{ℓm_ℓ}
ωℓ1,..,ωℓmℓ 并且我们使用 $x_{ℓj} 作为
x
(
ω
ℓ
j
)
x(ω_{ℓj})
x(ωℓj)的简写。 因此,单独的过程由下式给出
z
ℓ
(
s
)
=
∑
j
=
1
m
ℓ
k
ℓ
(
s
−
w
ℓ
j
x
ℓ
j
)
z_ℓ(s)=\sum_{j=1}^{m_ℓ}k_ℓ(s-w_{ℓj}x_{ℓj})
zℓ(s)=∑j=1mℓkℓ(s−wℓjxℓj)
多分辨率被捕获在随着
ℓ
ℓ
ℓ 增加而变窄的内核中。 因此,每个额外的
z
ℓ
(
s
)
z_ℓ(s)
zℓ(s) 都说明了额外的小尺度细节。 与基本公式一样,模型可以表示为混合模型
y
=
μ
1
+
∑
ℓ
=
1
p
K
ℓ
x
ℓ
+
ϵ
y=μ1+\sum_{ℓ=1}^pK^ℓx_ℓ+ ϵ
y=μ1+∑ℓ=1pKℓxℓ+ϵ
其中
x
ℓ
x_ℓ
xℓ 是
m
ℓ
m_ℓ
mℓ 向量
(
x
ℓ
1
,
.
.
.
x
ℓ
m
ℓ
)
T
(x_{ℓ1}, . . . x_{ℓm_ℓ} )^T
(xℓ1,...xℓmℓ)T ,
K
ℓ
K^ℓ
Kℓ 的元素由下式给出
K
i
j
ℓ
=
k
ℓ
(
s
i
−
ω
ℓ
j
)
x
ℓ
j
K^ℓ_{ij}=k_ℓ(s_i-ω_{ℓj})x_ℓj
Kijℓ=kℓ(si−ωℓj)xℓj
x
ℓ
x_ℓ
xℓ~
N
(
0
,
σ
x
ℓ
2
I
m
ℓ
)
N(0, σ^2_{x_ℓ}I_{m_ℓ})
N(0,σxℓ2Imℓ)
ϵ
ϵ
ϵ ~
N
(
0
,
σ
ϵ
2
I
n
)
N(0, σ^2_ϵI_n)
N(0,σϵ2In)
4.1 简单的一维示例
例如,图 7 显示了前一维示例中的 n = 30 个数据点。 空间位置
s
i
s_i
si在 1 到 10 之间变化,实际数据点是根据公式创建的:
y
(
s
i
)
=
s
i
n
(
2
∗
π
[
s
i
/
10
]
)
+
0.2
s
i
n
(
2
∗
π
[
s
i
/
2.5
]
)
+
e
i
y(s_i)=sin(2*π[s_i/10])+0.2sin(2*π[s_i/2.5])+e_i
y(si)=sin(2∗π[si/10])+0.2sin(2∗π[si/2.5])+ei
其中
e
i
s
e_is
eis是独立同分布
N
(
0
,
0.
1
2
)
N(0,0.1^2)
N(0,0.12),所以第一项给出了大规模的变化,而第二项给出了 4 倍频率的五分之一的变化。 多分辨率模型构造如下:
所以第一项给出了大规模的变化,而第二项给出了 4 倍频率的五分之一的变化。 多分辨率模型构造如下:
图 7 基本模型(实线)和多分辨率模型(虚线)下拟合值的模拟空间数据集。 在这里,多分辨率模型在数据中找到小尺度结构。
拟合过程显示为来自基本粗略模型的拟合。 请注意,多尺度公式从高频 s i n sin sin 项中提取信号。 下面给出了两种模型下的估计方差分量。 请注意,多分辨率模型以最高分辨率 ( ℓ = 3 ) (ℓ = 3) (ℓ=3) 获取高频 sin() 项,其中 k 3 ( s ) k_3(s) k3(s) 几乎匹配正弦波的半个周期。
4.2 一个二维示例
类似地,这种多分辨率模型可以应用于上一节的臭氧示例。 这里指定多分辨率模型具有
p
=
2
p = 2
p=2 个级别。 第一个粗分量
(
ℓ
=
1
)
(ℓ = 1)
(ℓ=1)正是原始臭氧模型的规范,因此第 3.2 节给出了对
x
1
(
s
)
,
k
1
(
s
)
{x_1(s), k_1(s)}
x1(s),k1(s)。 精细组件在两个空间坐标中将粗略规格细化为两倍。 因此,过程
x
2
(
s
)
x_2(s)
x2(s) 被限制在一组更密集的点
ω
2
,
1
,
.
.
,
ω
2
,
87
ω_{2,1}, . . , ω_{2,87}
ω2,1,..,ω2,87 显示在图 8 的左下角。除此之外,内核
k
2
(
s
)
k_2(s)
k2(s) 的尺度只有
k
1
(
s
)
k_1(s)
k1(s) 的一半,因此内核从其中心消失 4.5 度。
图 8. 应用于臭氧数据的基本(粗略)和多分辨率模型。 左上:潜在进程的
ω
1
j
(
s
)
ω_{1j}(s)
ω1j(s)的数据和位置; 右上:符合基本公式; 左下:细尺度过程对应的
ω
2
j
(
s
)
ω_{2j}(s)
ω2j(s)的数据和位置——粗略的
ω
1
j
(
s
)
ω_{1j}(s)
ω1j(s)如左上图所示; 右下:多分辨率公式拟合。
使用 REML 拟合数据集,结果总结如下。 预测值显示在图 8 的右侧框架中。与一维示例一样,多分辨率拟合与基本模型的不同之处在于更精细的尺度细节。 臭氧场的这种分解是有吸引力的,因为预计较小规模的波动是相当不可预测的。 然而,臭氧浓度的较大规模变化确实显示出每天的一些持久性。
5.建立时空模型
也许这些过程卷积模型的最大吸引力在于它们为开发新的空间和时空模型提供了一个框架,这些模型允许更真实的时空依赖性,同时保持一些分析的易处理性。 通常,可以通过首先定义一个简单的、可能是离散的空间和时间上的过程,然后用一个或多个内核对其进行平滑处理,从而构建一个空间和时间上的平滑过程,从而构建一个时空过程。
这种建设性的方法很有吸引力,因为得到的模型可以扩展以允许泛化,例如非平稳性、非高斯模型和不可分离的时空依赖结构。 参见 Wolpert 和 Ickstadt (1998)、Ickstadt 和 Wolpert (1999) 以及 Higdon 等人。 (1999) 用于一些纯粹的空间应用,Higdon (1998) 用于时空模型。 此外,可以以促进计算的方式构建模型——例如将底层过程限制在格上,以便可以使用快速傅里叶变换。
例如,不可分离的非高斯时空过程的一个简单但重要的示例可以构造如下:
- 在空间 s 和时间 t 上按照某个标记点过程分布随机变量。
- 根据在空间 s 和时间 t 上定义的一些内核来平滑这个过程。
图 9 显示了从这种过程生成实现的方案。在这个特定示例中,i.i.d. 指数随机变量与每个位置相关联,从而产生严格的正非高斯场。 请注意,这里的空间只是一维的,以使图形易于理解——也可以考虑 2 维或 3 维空间分量。
对于恒定的盛行风下的污染物浓度,这个相当简单的例子可能是令人满意的。 然而,它可能不能很好地解释不断变化的风型。 在这种情况下,最好让基础点过程随着时间的推移而演变,并且只在空间分量上平滑。 图 10 显示了在空间点阵上定义的潜在过程,该过程在幅度和空间位置上随时间缓慢变化。 这样一个过程很有吸引力,因为这个潜在过程的空间“转移”可能与该地区的气象数据有关。
图 9. 平滑在空间和时间上定义的标记点过程产生时空表面。
5.1 臭氧浓度的时空模型
利用上述思路,构建了连续 30 天臭氧浓度的时空模型。 我们让
t
∈
T
=
1
,
.
.
.
,
30
t ∈ T = {1, . . . , 30}
t∈T=1,...,30 索引时间(以天为单位)并定义潜过程的空间支持
W
=
ω
1
,
.
.
.
,
ω
27
W = {ω_1, . . . , ω_{27}}
W=ω1,...,ω27 与第 3.2 节的基本模型完全相同,因此
x
(
s
,
t
)
x(s, t)
x(s,t) 在集合
W
×
T
W × T
W×T 上非零。 而不是将随机变量
x
j
t
=
x
(
ω
j
,
t
)
x_{jt} = x(ω_j , t)
xjt=x(ωj,t) 定义为 i.i.d.,每个序列 {
x
j
t
x_{jt}
xjt},
t
=
1
,
.
.
.
,
30
t = 1, . . . , 30
t=1,...,30 被指定为遵循高斯随机游走。 通过考虑潜在过程
x
(
s
,
t
)
x(s, t)
x(s,t) 内的时间依赖性,通过在空间上平滑
x
(
s
,
t
)
x(s, t)
x(s,t) 获得诱导时空过程
z
(
s
,
t
)
=
∫
s
k
(
u
−
s
)
x
(
u
,
t
)
=
∑
j
k
(
ω
j
−
s
)
x
j
t
z(s,t)=\int_sk(u-s)x(u,t)=\sum_jk(ω_j-s)x_{jt}
z(s,t)=∫sk(u−s)x(u,t)=∑jk(ωj−s)xjt
数据由 30 天中每天记录的大约 500 次测量组成。 前 8 天的测量结果如图 11 所示。台站这么多,每天都有一小部分无法上报数据。 因此,每个时间步长
n
t
n_t
nt 记录的实际数据点数量会有所不同。 以潜在过程值
x
t
=
(
x
1
,
t
,
.
.
.
.
,
x
27
,
t
)
T
,
t
=
1
,
.
.
.
,
30
x_t = (x_{1,t}, . . . . , x_{{27},t})^T , t = 1, . . . , 30
xt=(x1,t,....,x27,t)T,t=1,...,30, 数据
y
t
=
(
y
1
t
,
.
.
.
,
y
n
t
t
)
T
y_t = (y_{1t}, . . . , y_{{n_t}t})^T
yt=(y1t,...,yntt)T 的模型,记录在站点
s
1
t
,
.
.
.
,
s
n
t
t
s_{1t}, . . . , s_{{n_t}t}
s1t,...,sntt 可以表示为两个演化方程
y
t
=
K
t
x
t
+
ϵ
t
y_t=K^tx_t+ϵ_t
yt=Ktxt+ϵt(6)
x
t
=
x
t
−
1
+
v
t
x_t=x_{t-1}+v_t
xt=xt−1+vt(7)
其中,
K
t
K^t
Kt是
n
t
×
27
n_t×27
nt×27的矩阵,由下式给出
K
i
j
t
=
k
(
s
i
t
−
ω
j
)
x
j
t
,
t
=
1
,
.
.
.
30
,
K^t_{ij}=k(s_{it}-ω_j)x_{jt},t=1,...30,
Kijt=k(sit−ωj)xjt,t=1,...30,
ϵ
t
ϵ_t
ϵt ~
(
i
i
d
)
(iid)
(iid)
N
(
0
,
σ
ϵ
2
)
,
t
=
1
,
.
.
.
30
,
N(0,σ^2_ϵ),t=1,...30,
N(0,σϵ2),t=1,...30,
v
t
v_t
vt~
(
i
i
d
)
(iid)
(iid)
N
(
0
,
σ
v
2
)
,
t
=
1
,
.
.
.
30
,
N(0,σ^2_v),t=1,...30,
N(0,σv2),t=1,...30,
x
1
x_1
x1~
N
(
0
,
,
σ
x
2
I
27
)
N(0,,σ^2_xI_{27})
N(0,,σx2I27)
该模型很容易适应 West 和 Harrison (1997) 的动态线性模型 (DLM) 机制。 其他替代方法是通过 MCMC 或基于 REML 的方法进行完全贝叶斯分析。 参见 Stroud 等人。 (1999) 用于非常相似的基于 DLM 的建模方法。
对于此示例,为精度
1
/
σ
ϵ
2
1/σ^2_ϵ
1/σϵ2 、
1
/
σ
v
2
1/σ^2_v
1/σv2 和
1
/
σ
x
2
1/σ^2_x
1/σx2 指定了平坦但适当的 Gamma 先验。 然后使用 MCMC 进行估计。拟合值的后验均值估计值如图 12 所示。此外,图 13 显示了
x
(
s
,
t
)
x(s,t)
x(s,t) 在三个内部空间支持点处作为时间函数的后验均值。
6 讨论
本文重点关注构建空间和时空模型的工具,而不是关注实际的数据分析。 希望让读者了解这些模型如何用于应用程序。
该方法在精神上类似于使用经验正交函数来降低空间场的维数(例如,参见 von Storch 和 Zwiers,1999 年),但过程卷积方法提供了一些优势。 像小波一样,过程卷积允许对空间场进行局部控制。 插值是“内置的”来处理卷积; 应该如何插入数据的 EOF 表示不太明显。 最后,该方法可以通过标准的高斯过程理论来理解。 因此,关于构建高斯过程模型的大量文献也适用于过程卷积模型。
最后,我注意到关于臭氧数据建模的研究正在进行中。 一些令人感兴趣的问题是:
1.日常臭氧场非平稳性的本质是什么? 如果依赖性确实随空间位置而变化,那么这种非平稳性的性质是否每天都相似? 尼奇卡等人。 (1999) 给出了一种基于小波的方法,用于从随时间的重复观察中估计非平稳性。
2.在多分辨率模型中,有多少分辨率级别是合适的? 据推测,这个问题的答案取决于实际的推理问题,以及要建模的空间过程。 此外,数据可能只提供有限的分辨率。 任何有关较小规模的信息都必须从其他来源获得并内置到先前的模型中。 如此精细的尺度信息如何获得?
3.将多分辨率模型与时空模型结合起来可能是有利的。 空间模型的粗分辨率组件可能显示出时间依赖性,而更精细的分辨率组件则没有。 应该如何确定应该将哪些决议纳入模型的时间组件?
4.通过调节几个估计参数,可以得到封闭形式的后验分布。 这种后验形式将适用于强大的基于效用的设计方法(Müler,1999)。 了解这种“条件”后验与更合适的边缘后验有何不同将是很有趣的。
我们希望这一系列研究能够在进一步了解臭氧和其他污染物的性质和可变性方面发挥作用。