【摘 要】 本文简要介绍了贝叶斯分层建模方法的概念、优势和局限性。
【原 文】 N. Cressie, chapter 2, Statistics for spatio-temporal data. 2011.
核心内容快速浏览
(1)贝叶斯全概率公式
贝叶斯全概率公式允许将随机变量的联合分布分解为一系列条件分布:
[A,B,C]=[A∣B,C][B∣C][C]
其中 “[⋅]” 用于表示概率分布;例如,[A,B,C] 是随机变量 A、B 和 C 的联合分布,而 [A∣B,C] 是给定 B 和 C 时 A 的条件分布。
(2)Berlinear 的贝叶斯分层模型 (BHM) 范式
Mark Berliner (Berliner,1996)是最早使用贝叶斯全概率公式分解来为复杂过程建模的人。也就是说,联合分布 [data,process,parameters],可以自顶向下分解为三个层次:
数据模型:在假设下层的 “真实” 过程(有时称为隐过程或潜在过程)已知,并给定一些参数的情况下,描述数据及其不确定性的一个概率模型。
过程模型:在给定一些参数的情况下,描述隐过程(及其不确定性)的一个概率模型。在此级别中,模型不需要考虑测量不确定性(即偶然不确定性)。过程模型可以基于科学理论或经验知识,并且这些知识通常是物理性的或机械性的。
参数模型:对数据模型和过程模型中的参数的不确定性进行建模。
从上到下,贝叶斯分层模型的级别是:
- 数据模型:[data | process, parameters]
- 过程模型:[process|parameters] parameters]
- 参数模型: [parameters]
重要的是,上述每一个级别都可以有子级别,只要能够给出条件概率模型。
事实上,我们最终感兴趣是后验分布:[process, parameters∣data],而最方便的地方在于,我们关心的后验分布与上面给出的 贝叶斯分层模型各层乘积(未归一化的似然)成正比:
[process,parameters | data]∝××[data∣process,parameters][process∣parameters][parameters]
其中“ ∝ ”表示 “与某某成正比”。事实上,根据贝叶斯定理,将上式右侧除以归一化常数 [data],就可以等于左侧。但众所周知,[data] 很难有闭式解,而且给定数据的情况下是个常数,所以很多时候使用正比符号来做相对计算和比较。当确实需要推断得出后验分布时,可以采用一些数值方法来(例如 MCMC 方法)实现,贝叶斯模型最近几十年的快速发展,正是得益于这些数值方法的提出。
(3)经验分层模型 (EHM) 范式
经验分层模型 (EHM) 仅使用前两个级别,预测分布为:
[process | data,parameters]∝×[data|process, parameters][process| parameters]
贝叶斯分层模型中作为随机项的 parameters 在 EHM 中被替换为根据经验数据或知识获得的参数估计值(即确切的已知值,这种情况在某些应用领域也很常见)。对于 EHM,不需要对参数做推断,但依然需要考虑预测分布的推断问题,同样出于无闭式解的原因,需要通过 MCMC 等方法对预测分布进行数值估计。
1 术语解释
(1)统计推断:
数据、过程或参数的不确定性意味着结论存在不确定性。统计学家们将这种在存在不确定性的情况下得出的结论称为统计推断(或推断);在本书中,推断要么是对固定但未知参数的估计(Estimation),要么是对未知随机量的预测(Prediction)。请注意,“Forecasting” 指对未来做出结论,可以被视为 “Prediction” 的一个特例。
(2)正态分布和高斯分布
两者是同义词,但本文更倾向于后者。
使用表达式 Z∼Gau(μ,σ2) 来表示服从高斯分布的随机变量 Z,其均值为 μ、方差为 σ2;也可以等价地表示为 Z∼N(μ,σ2)。
随机向量 Z≡(Z1,…,Zm)T 是一个 m 维的列向量,其中符号 “ T ” 表示转置。则 Z∼Gau(μ,Σ) 表示均值向量为 μ、协方差矩阵为 Σ 的 m 维高斯分布。
协方差矩阵 Σ 有时也被称为 “方差矩阵” 或 “方差-协方差矩阵”,是一个 m×m 的对称正定矩阵(有时允许非负定),其第 (i,j) 个元素的值为 Cov(Zi,Zj);i,j=1,…,m。
(3)概率分布的表示
可以令 [A] 代表随机量 A 的概率分布(如果不习惯,也可以等价地表示为 p(A))。例如:
单随机变量 Z∼Gau(μ,σ2) 可以等价地表示为:
[Z]=(2πσ2)−1/2exp[−(1/2)(z−μ)2/σ2]:z∈R
多变量(随机向量) Z∼Gau(μ,Σ) ( Z=(Z1,…,Zm)T )可以表示为:
[Z]=(2π)−m/2∣Σ∣−1/2exp[−(1/2)(Z−μ)TΣ−1(Z−μ)]:z∈Rm
(4)概率分布的性状描述
如果随机变量 X 的概率分布可以被参数化为函数形式 f,则可以使用术语 “矩” 来描述分布的特征:
原点矩:
零阶矩
为总概率 1;
一阶矩
为均值 (Mean),反映了分布的中心位置,定义为 Mean(X)=μ=E[X]=∫−∞+∞xf(x)dx。一阶矩是基于原点的矩类型,反映了相对于原点的位置特征。
中心矩:
在一阶矩基础上,通常采用以均值为中心的中心矩方式来定义其他矩,以便更好地描述分布的形状信息。
一阶中心距
显然为 0;
二阶中心矩
为方差 (Variance),反映了分布的离散度(或反之,聚集程度),定义为 Var(X)=σ2=E[(x−E[X])2]=∫−∞+∞(x−μ)2f(x)dx;
在二阶中心矩基础上,可以对更高阶的矩做归一化处理,得到高阶的标准化矩:
三阶标准化距
为偏度 (Skewness),反映了分布的对称程度,定义为 Skew(X)=E[(σX−E[X])3]=∫−∞+∞(σx−μ))3f(x)dx 。偏度 =0 时,为对称分布;偏度 >0 时为右偏分布;偏度 <0 时为左偏分布。
四阶标准化距
为峰度 (Kurtosis),反映了分布中峰值的尖度,定义为 Kurt(X)=E[(σX−E[X])4]=∫−∞+∞(σx−μ))4f(x)dx
多变量的混合矩
在单变量矩以及中心矩基础上,可以定义两个随机变量 X 和 Y 之间的 k+p 阶混合中心距 E[(X−E[X])k(Y−E[Y])p] 。
常用的混合中心矩为(1+1)阶混合中心距,即协方差 Cov(X,Y)=E[(X−E[X])(Y−E[Y)]
根据矩的定义,当概率密度函数 f 比较简单时(如常见的高斯等分布),能够得到矩的解析形式。但现实大部分场景中,随机变量的分布是复杂的,通常无法得到解析解,只能通过样本来获得近似的数值解。有了对矩的认识,我们就可以从样本数据中计算矩,以使两者的充分统计量的期望能够匹配。
(5)随机变量的函数
如果 g(A) 是由函数 g(⋅) 明确定义的随机量,那么其期望值 E[g(A)],在连续情况下,可以等价地写为 E[g(A)]=∫Ag(A)[A]dA;在离散情况下可以等价地写为 E[g(A)]=Σg(A)[A]。离散形式为函数 g(A) 的蒙特卡洛数值积分提供了支撑。
(6)条件分布
令 [A∣B] 表示在指定随机变量 B 的特定值的条件下,随机变量 A 的条件分布。简称为给定 B 时 A 的条件分布。例如:
表达式 Z∣Y=y∼Gau(y,σ2I) 可以等效地写成:
[Z∣Y=y]=(2πσ2)−m/2exp[−(1/2)(z−y)T(z−y)]:z∈Rm
(7)时空索引与时空模型
令我们的感兴趣空间域为 Ds⊂Rd,是 d 维欧几里得空间的一个子集,令我们感兴趣的时间域为 Dt⊂R1。空间索引 s(一个 d 维向量)和时间索引 t (实数)可以在其各自域 Ds 和 Dt 上连续或离散变化。
本文只关注时空随机过程。当时间 t 连续变化时,我们将随机过程写为:
{Y(s;t):s∈Ds,t∈Dt};
当 t 为离散变化时,我们使用下面的形式以示区别。
{Y(s):s∈Ds,t∈Dt}
2 分层建模方法
2.1 分层模型
有一种被称为分层(统计)建模的通用方法可以用于表达不同来源的不确定性。本文主题是时空数据统计,因此感兴趣的对象主要与 “随机变量”、“随机向量” 或 “随机集” 的时空随机过程有关。
分层模型(缩写为 HM)中指定的条件概率分布通常取决于未知的参数:
- 贝叶斯分层模型:如果为了概率性地表达参数上的不确定性,在分层模型中包含了参数模型,则此时的分层模型被称为贝叶斯分层模型(BHM)。
- 经验分层模型:如果假设参数是固定的,只是还不知道具体值,需要使用数据来估计它们;然后再代入数据模型和流程模型中,就好像参数是已知的一样,则此时的分层模型被称为经验分层模型(EHM)。
本文中,术语 “先验分布” 与 “参数模型” 同义。
HM 中三个通用量是 Z、Y 和 θ,我们通常认为它们都是随机变量。其中将 Z 被视为数据,将 Y 被视为我们希望预测的(隐)过程,将 θ 被视为未知参数。Z、Y 和 θ 在现实案例中可能很复杂,例如:对给定一周内区域空气质量的空间统计制图,Z 可能是 100 维的向量,Y 可能是 1000 维的向量,而 θ 可能是五维。我们希望根据 Z 对 Y 和 θ 进行推断。也就是说,在 贝叶斯分层模型中,我们希望预测 Y 和 β 的分布,而在 EHM 中,我们希望预测 Y 的分布并估计出参数 β 的值。
2.2 贝叶斯分层模型
贝叶斯分层建模(Bayesian Hierarchical Modeling,BHM)的基本表示法将模型分成三个级别(Berliner,1996):
data model process model parameters model :[Z∣Y,θ]:[Y∣θ:[θ]
有时我们会用 [Z∣Y,θD] 和 [Y∣θP] 来强调和区分数据模型参数 θD 和过程模型参数 θP,则此时 θ={θD,θP},参数模型相应地表示为 [θD,θP]。则三者的联合概率分布为:
[Z,Y,θ]=[Z,Y∣θ][θ]=[Z∣Y,θ][Y∣θ][θ]
也就是说,联合分布是数据模型、过程模型和参数模型的乘积。联合分布的一种特殊情况是当 θ=θ0 已知时, [θ] 将其所有概率质量集中在点 θ0 处,此时对应于经验分层模型。
根据贝叶斯定理,当数据已知时,我们的目标是推断出过程模型 Y 和参数模型 θ 的后验分布。即:
[Y,θ∣Z]=∫∫[Z,Y,θ]dYdθ[Z∣Y,θ][Y,θ]=[Z][Z∣Y,θ][Y∣θ][θ]
在贝叶斯决策理论的框架内,贝叶斯分层模型中对 Y 和 θ 的所有推断都依赖于上述分布。
3 分层建模方法的优势
当依赖关系变得复杂时,HM 的真正威力才会变得显而易见。因为如有必要,人们可以进一步分解每个分量中的分布,并且可以通过建模假设对其进行简化。在下文中,我们将介绍并讨论分层模型的优势和局限性。
3.1 数据模型
假设我们对一个过程 Y 感兴趣,我们有几个不同的数据集 Z1、Z2、Z3,这些数据集都是对过程 Y 的含不确定性的测量值,也许它们是在不同空间或时间尺度上测量的。此时通常可以做出以下数据建模假设:
[Z1,Z2,Z3∣Y,θD]=[Z1∣Y,θD,1][Z2∣Y,θD,2][Z3∣Y,θD,3]
其中数据模型参数 θD 由 {θD,1,θD,2,θD,3} 给出。也就是说,我们可以假设不同数据集之间相互独立,这取决于真实的过程。
尽管科学必须证明上述假设的合理性,但通常它是成立的,并且为综合考虑各类型观测提供了一种非常方便的手段。每个子分布中的参数,可以拟合不同分辨率和对齐上的变化,以及不同的测量误差特性(例如,Gelfand、Zhu 和 Carlin,2001;Wikle 和 Berliner,2005)。
也可以对不同过程(或参数)的不同数据集进行组合。例如,假设有两个数据集 Z1 和 Z2,分别对应于过程 Y1 和 Y2,其中 Y=(Y1,Y2)。我们通常可以使用条件独立假设,如下所示:
[Z1,Z2∣Y,θD]=[Z1∣Y,θD,1][Z2∣Y,θD,2]
这里的数据模型参数由 θD={θD,1,θD,2} 给出。
3.2 过程模型
过程模型的分布也是可以分解的。例如,假设 Y 由两个子过程 Y1 和 Y2 组成,我们通常可以使用条件概率建模:
[Y1,Y2∣θP]=[Y1∣Y2,θP][Y2∣θP]
其中 θP 是过程模型的参数。也就是说,通常可以通过使用条件概率来简化过程的两个组分的联合交互。
最具代表性且众所周知的例子是具有马尔可夫性的时间序列模型,具体来说,考虑时间序列,Y1,Y2,…,YT–1,YT,通常很难指定所有这些随机变量的联合分布,但如果考虑对其进行条件概率分解:
[YT,…,Y1]=[YT∣YT−1,…,Y1][YT−1∣YT−2,…,Y1],…,[Y2∣Y1][Y1]
并对其做简化,则有可能构造出有价值的联合分布。
例如,使用一阶马尔可夫假设,即 t 时刻的过程仅取决于 t−1 时刻的过程,而和之前的时刻无关:
[Yt∣Yt−1,Yt−2,…,Y2]=[Yt∣Yt−1]
此时联合分布可以简化为:
[YT,…,Y1]=[YT∣YT−1[YT−1∣YT−2],…,[Y2∣Y1][Y1]
请注意,右侧现在仅取决于单变量和双变量的概率分布。上述性质在卡尔曼滤波的更新机制中起到决定性作用。Meinhold 和 Singpurwalla (1983) 证明其更新机制是对隐时间序列 {Yt:t=1,2,...} 进行时间分层建模的优化预测结果。当 Yt 是一个随机空间过程时,可以构建时空分层模型,从而产生时空卡尔曼滤波器(Huang 和 Cressie,1996;Wikle 和 Cressie,1999)。
这种过程模型的分解能力也同样很具有吸引力,因为它们利用了科学上合理的建模假设。非常复杂的联合概率分布可以通过相对简单的条件概率定义来建模。在过程模型中,确定性模型可以用随机成分重新指定以解释科学中存在的不确定性。例如,Wikle (2003b) 使用 “响应-扩散偏微分方程” 来驱动一个入侵物种的过程模型。
过程模型也可以进一步分解或从基于科学的建模假设中简化。例如,Wikle (2003b)以栖息地协变量和空间随机场为条件,指定了 “空间显式扩散系数的参数分布”。 人们也可以对后面的过程进行建模(以参数为条件),以将子层级合并到过程模型中。分层模型可以有很多层次,只要有科学见解或数据作为分解的基础。当不再有科学见解时,贝叶斯分层模型通常会采用 “非信息性” 先验作为参数模型的一部分。
3.3 参数的推断
即便大家同意分层模型的合理性,也仍然存在如何对其进行推断的问题。当然,我们可以只在前两个层次([Z∣Y,θD] 和 [Y∣θP] )上对不确定性进行建模,参数 θ={θD,θP} 可以认为是固定但未知的。这就是经验分层建模方法。 经典的线性混合模型可以被认为是经验分层模型(例如,Demidenko,2004)。空间预测(克里金法)与卡尔曼滤波器等顺序时间序列方法(例如,Meinhold 和 Singpurwalla,1983 年)一样适用于该框架(Cressie,1988)。
根据数据模型和过程模型的复杂性,通常可以使用经典统计估计方法来获得参数 θ 的估计。在经验分层模型上下文中,常用的估计方法包括最大似然估计、期望最大化(EM)算法、条件似然法、伪似然法以及估计方程(Demidenko,2004)。虽然经验分层模型没有明确认识到参数估计的不确定性,但这种不确定性通常可以通过重采样和自举程序来解释(Stoffer 和 Wall,1991;Efron 和 Tibshirani,1993;Lahiri,2003)。 Cressie (1993, pp. 492–497) 讨论了空间和时间上下文中的自举(包括参数自举)。在简单的情况下,可以使用统计性的扰动参数来近似地解释不确定性(Rao,2003,第 6.2 节)。
在复杂的分层模型中进行推断的一种相干方法是通过添加参数模型 [θD,θP] 来解释参数值的不确定性,这就是贝叶斯分层模型。回想一下,在给定数据的情况下,推断是基于过程和参数的后验分布:
[Y,θD,θP∣Z]∝[Z∣Y,θD][Y∣θP][θD,θP]
在实际工作中,很少有可能获得该后验分布的解析表达式。马尔可夫链蒙特卡罗 (MCMC) 计算方法可用于模拟一般贝叶斯分层模型的新发现(Gelfand 和 Smith,1990 年),彻底改变了统计科学,并且使更复杂的建模场景成为可能。
4 局限性
在过程模型中确定因果关系的问题可以用条件概率来表示。如果 Y1 是一种可以通过物理/化学/生物/经济等机制直接影响 Y2 的现象,并且 [Y2∣Y1] 随着 Y1 的变化而变化,则 Y1 是 Y2 原因的候选者。
但即便是最好的理论,也可能会遗漏一个重要因素(例如,F),这可能会削弱关系,甚至将正依赖关系逆转成了负依赖关系。我们将此称为 辛普森悖论
,而 F 称为 “潜伏变量”。
对于一个简单过程 Y=(Y1,Y2,F),过程模型可以写成
[Y]=[Y1,Y2∣F][F]=[Y2∣Y1,F][Y1∣F][F]
由此,对于问题 “为什么 Y2 会有如此行为?” ,可以有这样的答案: “ 因为 Y1 导致了这种行为,而这种依赖的类型受因素 F 控制”。如果错误地聚焦在 [Y2∣Y1] 而不是 [Y2∣Y1,F] 上,我们可能会推断出不正确的依赖关系。当 F 代表 “空间聚合水平” 时,此类虚假推断被称为 生态学谬误(Robinson,1950);对此的空间统计解释在 Cressie (1996) 中可以找到,被称为 支撑的变化。
更一般地说,许多因素具有空间和时间可变性;因此,如果过程模型无法解释这种可变性,空间和时间就完全可以扮演 F 的角色。换言之,对时空可变性建模以及良好的实验设计,可以让我们更接近科学的真相 – 因果关系。
5 其他
建立(贝叶斯)分层统计模型然后进行有效推断存在许多挑战。对贝叶斯方法的历史批评是:它们需要对参数的先验进行 “主观” 指定。当然,经典的边缘概率模型中对似然的指定也存在主观性。事实上,更广泛的观点是:所有模型组件的指定都涉及主观性,这就包括数据模型、过程模型和参数模型。但 “主观” 的含义并不总是很清晰。使用确定性关系来驱动随机模型可能是 “主观的”,例如热带风(Wikle 等,2001),但这种模型所依据的科学来自牛顿定律运动!因此,我们认为将判断统计模型的概率分布的主观/客观性当作一个问题,实际上没有任何意义。精力应当更多放在对模型选择推断的敏感性上,还有就是该选择在科学上是否有意义。
鉴于在开发模型时建模者会面临大量信息,而分层模型的条件概率框架可以用来识别这些信息。例如 I 是条件作用的一部分。对于贝叶斯分层模型,我们有:
[Y,θD,θP∣Z,I]∝[Z∣Y,θD,I][Y∣θP,I][θD,θP∣I]
此范式中面临的一个主要挑战是:如何在可能的范围内了解信息 I 的重要性?通常情况下,一组研究人员有一个集体 “I”,并且会比任何个人的 “I” 都量化的更合适。