《A New Look at the Statistical Model Identification》
统计模型识别的新视角
转载自 IEEE Transactions on Automatic Control,第 19 卷,第 716-723 页,1974。
作者: 赤池弘嗣,IEEE 会员
【摘要】
本文简要回顾了时间序列分析中统计假设检验的发展历史,并指出将假设检验过程简单等同于统计模型识别的方法是不充分的。文中回顾了经典的最大似然估计过程,并引入了一种为统计识别目的而设计的新估计——最小信息理论准则(AIC)估计(MAICE)。当存在多个竞争模型时,MAICE 被定义为:选择使下式最小的模型及其参数的最大似然估计值
AIC
=
(
−
2
)
log
(最大似然)
+
2
(
模型中独立调节参数的个数
)
\operatorname{AIC} = (-2) \log \text{(最大似然)} + 2 \, (\text{模型中独立调节参数的个数})
AIC=(−2)log(最大似然)+2(模型中独立调节参数的个数)
MAICE 提供了一种多功能的统计模型识别程序,其不受传统假设检验方法中固有的歧义影响。文中以一些数值例子演示了 MAICE 在时间序列分析中的实际应用效用。
1.引言
尽管统计概念和模型最近已在工程和科学的几乎每个领域中得到发展,但基于有限数量观测所提供的信息构造一个恰当模型的困难似乎尚未被完全认识。毫无疑问,统计模型构造或识别这一课题在很大程度上依赖于对被观察对象理论分析的结果。然而,必须意识到理论结果与实际识别程序之间通常存在很大的差距。一个典型的例子是线性系统最小实现理论的结果与基于有限记录识别随机过程的马尔可夫表示之间的差距。线性系统的最小实现通常是通过对某个 Hankel 矩阵的行或列的秩或依赖关系进行分析而定义的 [1]。在实际情况下,即使该 Hankel 矩阵在理论上是给定的,由于舍入误差,矩阵总是表现为满秩;如果该矩阵是从真实物体的观测记录中获得的,那么矩阵各元的抽样变异性远大于舍入误差,而且系统总是无限维的。因此,可以看出,统计识别问题本质上涉及对真实结构的近似,这是一种人类智力活动的基本艺术。
正如 Lehman [2, p. viii] 指出,假设检验程序传统上应用于实际上需要多个决策程序的情况。如果把统计识别程序看作一种决策程序,那么最根本的问题在于如何适当地选择损失函数。在 Neyman-Pearson 统计假设检验理论中,只考虑接受正确与错误假设的概率来定义决策所造成的损失。在实际情况中,所假定的原假设仅是近似,并且几乎总是与实际情况不同。因此,在检验理论中损失函数的选择使得其实际应用在逻辑上具有矛盾性。认识到这一点——即假设检验程序并未充分被表述为一种近似程序——对发展具有实际应用价值的识别程序非常重要。
通过对极具实用性和成功性的最大似然方法的分析,可以获得对识别问题的新视角。最大似然估计在某些正则条件下是渐近有效的,这表明似然函数往往是一种对参数在真实值附近的微小变化最为敏感的量。这一观察提示我们可以采用
S ( g ; f ( ⋅ ∣ θ ) ) = ∫ g ( x ) log f ( x ∣ θ ) d x S(g ; f(\cdot \mid \theta))=\int g(x) \log f(x \mid \theta) \, d x S(g;f(⋅∣θ))=∫g(x)logf(x∣θ)dx
作为衡量由概率密度函数
f
(
x
∣
θ
)
f(x \mid \theta)
f(x∣θ) 定义的概率结构的模型与由密度函数
g
(
x
)
g(x)
g(x) 定义的结构之间“拟合度”的准则。与经典最大似然估计过程中只考虑单一族的
f
(
x
∣
θ
)
f(x \mid \theta)
f(x∣θ) 的假定相反,在识别实际问题中通常会考虑几个替代模型或族,这些模型可以由不同形式的密度函数定义,和/或由相同形式但对参数向量
θ
\theta
θ 限制不同的情况定义。对最大似然估计(MLE)的详细分析自然而然引出了对于这种多模型情形下有用的新估计的定义。这一新估计称为最小信息理论准则(AIC)估计(MAICE),其中 AIC 代表作者最近所提出的信息理论准则 [3],它是对模型拟合优度的一种估计。MAICE 定义为:选择在各个模型中使 AIC(其定义为
AIC
=
(
−
2
)
log
(最大似然)
+
2
(
独立调节参数个数
)
\operatorname{AIC} = (-2) \log \text{(最大似然)} + 2 \, (\text{独立调节参数个数})
AIC=(−2)log(最大似然)+2(独立调节参数个数)
取得最小值的模型及其对应参数值。通过引入 MAICE,将统计识别问题明确表述为一个估计问题,从而完全消除了传统假设检验程序中为确定显著性水平而所需的主观判断。本文主要目的在于给出 MAICE 的明确定义,并通过与基于估计和假设检验的传统识别程序作比较,讨论其特性。虽然 MAICE 提供了一种可以在所有统计模型构建领域中使用的多用途识别方法,但其在时间序列分析中的实际效用尤为显著。文中给出一些数值例子,以显示 MAICE 如何为时间序列分析问题提供客观定义的答案,而相比之下传统假设检验方法只能给出主观且经常是不确定的答案。
2.时间序列分析中的假设检验
时间序列检验程序的研究始于对单个序列相关系数等于 0 这一简单假设的检验。显然这种检验的实用性过于有限,不能使其成为普遍有用的模型识别程序。1947 年,Quenouille [4] 引入了一种对自回归(AR)模型拟合优度的检验。Wold [5] 将 Quenouille 检验的思想推广到对移动平均(MA)模型拟合优度的检验。此后,对这些检验程序进行了若干改进和推广 [6]-[9],但是对时间序列假设检验领域做出最重要贡献的是 Whittle [10], [11],他系统地将 Neyman-Pearson 似然比检验程序应用于时间序列的情况。
时间序列中一个非常基本的检验是白噪声检验。在许多识别问题中,在拟合模型之后要求残差序列是白噪声,以证明模型的充分性,而白噪声检验在实践中被广泛使用 [12]-[15]。对于白噪声的检验,周期图分析提供了一种通用解法。Hannan [16] 对包括基于周期图检验在内的经典假设检验程序做了很好的阐述。
AR 或 MA 模型的拟合本质上是一个多重决策问题,而非单纯假设检验问题。Anderson [17] 明确讨论了 Gaussian AR 过程阶数的确定问题,将其作为一个多重决策程序。该程序采用了从最高阶开始依次向下至最低阶的模型序列检验形式。实际应用该程序时,必须为每个阶数的模型指定检验的显著性水平。虽然该程序设计上满足某些明确定义的最优性条件,但阶数判定问题的根本困难仍然在于如何选择各阶模型检验的显著性水平。同时,该决策程序的损失函数是由错误决策的概率定义的,因此该程序同样不免存在以下逻辑矛盾,即在实际应用中真实结构的阶数总是无限的。只有将问题明确重构为用模型对真实结构近似的问题,才能避免这一困难。
3.直接控制模型误差的方法
在非时间序列回归分析领域中,Mallows 为变量选择引入了一统计量 C p C_{p} Cp [18]。 C p C_{p} Cp 定义为
C p = ( σ ^ 2 ) − 1 (残差平方和) − N + 2 p C_{p}=\left(\hat{\sigma}^{2}\right)^{-1} \text{(残差平方和)}-N+2p Cp=(σ^2)−1(残差平方和)−N+2p
其中 σ ^ 2 \hat{\sigma}^{2} σ^2 是对真实残差未知方差 σ 2 \sigma^{2} σ2 的一个适当估计, N N N 为观测数,而 p p p 为回归中变量的个数。如果拟合的模型是正确的,则 C p C_{p} Cp 的期望大致为 p p p,否则将更大。 C p C_{p} Cp 是在使用估计的回归系数进行预测时,对预测平方误差(以 σ 2 \sigma^{2} σ2 标度)的期望和的一个估计,并被明确地解释为衡量所采用模型充分性的标准。由于 C p C_{p} Cp 在定义中对 σ ^ 2 \hat{\sigma}^{2} σ^2 的选择需要一些主观判断,因此引起了关注实际数据回归分析者的重视(参见文献 [18])。不幸的是, C p C_{p} Cp 的定义中对 σ ^ 2 \hat{\sigma}^{2} σ^2 的选择需要一定的主观判断。
几乎在 C p C_{p} Cp 被引入的同时,Davisson [19] 分析了当使用预测器的估计系数进行预测时平稳 Gaussian 过程的均方预测误差,并讨论了自适应平滑滤波器的均方误差问题 [20]。观测的时间序列 x i x_{i} xi 被视为信号 s i s_{i} si 与加性白噪声 n i n_{i} ni 的叠加。滤波输出 s ^ i \hat{s}_{i} s^i 定义为
s ^ i = ∑ j = − M L β j x i + j , ( i = 1 , 2 , ⋯ , N ) \hat{s}_{i}=\sum_{j=-M}^{L} \beta_{j} x_{i+j}, \quad (i=1,2,\cdots, N) s^i=j=−M∑Lβjxi+j,(i=1,2,⋯,N)
其中 β j \beta_{j} βj 是由样本 x i ( i = 1 , 2 , ⋯ , N ) x_{i} \, (i=1,2,\cdots, N) xi(i=1,2,⋯,N) 确定的。问题在于如何定义 L L L 和 M M M 以使得 N N N 个样本上的均方平滑误差
E [ 1 N ∑ i = 1 N ( s i − s ^ i ) 2 ] E\left[\frac{1}{N}\sum_{i=1}^{N}\left(s_i-\hat{s}_i\right)^2\right] E[N1i=1∑N(si−s^i)2]
达到最小化。在对 s i s_i si 和 n i n_i ni 作出适当假定之后,Davisson [20] 得到了这一误差的一个估计,其定义为
σ ^ N 2 [ M , L ] = s 2 + 2 c ^ ( M + L + 1 ) / N \hat{\sigma}_{N}^{2}[M, L]=s^2+2\hat{c}(M+L+1)/N σ^N2[M,L]=s2+2c^(M+L+1)/N
其中 s 2 s^2 s2 是对误差方差的估计,而 c ^ \hat{c} c^ 是在 ( M + L ) / N (M+L)/N (M+L)/N 的“较大”值处 s 2 s^2 s2 关于 ( M + L ) / N (M+L)/N (M+L)/N 曲线的斜率。该结果与 Mallows 的 C p C_{p} Cp 密切对应,并提示了这种统计量在用于预测的模型识别领域中的重要性。和 Mallows 的 C p C_{p} Cp 中对 σ ^ 2 \hat{\sigma}^{2} σ^2 的选择类似,当前统计量中对 c ^ \hat{c} c^ 的选择在实际应用中同样成为了一个难题。
在 1969 年,作者在不知与上述两程序密切关系的情况下,引入了一种一元 AR 模型拟合程序,其模型定义为
y i = a 1 y i − 1 + ⋯ + a p y i − p + x i , y_{i}=a_{1} y_{i-1}+\cdots+a_{p} y_{i-p}+x_{i}, yi=a1yi−1+⋯+apyi−p+xi,
其中 x i x_i xi 为白噪声 [21]。在这一程序中,通过控制使用最小二乘估计得到的系数进行一步预测的均方误差来实现目标。该均方误差被称为最终预测误差(FPE),当给定数据 y i ( i = 1 , 2 , ⋯ , N ) y_i \, (i=1,2,\cdots,N) yi(i=1,2,⋯,N) 时,其估计定义为
FPE ( p ) = { ( N + p ) / ( N − p ) } ⋅ ( C ^ 0 − a ^ p 1 C ^ 1 − ⋯ − a ^ p p C ^ p ) \begin{aligned} \operatorname{FPE}(p)&=\{(N+p)/(N-p)\} \\ &\quad\cdot\left(\hat{C}_{0}-\hat{a}_{p1}\hat{C}_{1}-\cdots-\hat{a}_{pp}\hat{C}_{p}\right) \end{aligned} FPE(p)={(N+p)/(N−p)}⋅(C^0−a^p1C^1−⋯−a^ppC^p)
其中假定 y i y_i yi 的均值为 0 0 0,且
C ^ l = ( 1 / N ) ∑ i = 1 N − l y i + l y i , \hat{C}_{l}=(1/N)\sum_{i=1}^{N-l} y_{i+l} y_{i}, C^l=(1/N)i=1∑N−lyi+lyi,
a ^ p i \hat{a}_{pi} a^pi 由由 C ^ i \hat{C}_{i} C^i 构成的 Yule-Walker 方程求解而得。通过将 p p p 从 0 依次扫描至某一上限 L L L,被识别出的模型是使 FPE ( p ) \operatorname{FPE}(p) FPE(p) (对于 p = 0 , 1 , ⋯ , L p=0,1,\cdots, L p=0,1,⋯,L)取得最小值的那个 p p p 及相应的 a ^ p i \hat{a}_{pi} a^pi。在这一程序中,除了上限 L L L 的确定需要主观判断以外,其它部分均没有留下主观成分。该程序的特性经过进一步分析 [22] 后,在实际数据中取得了非常出色的效果 [23], [24]。Gersch 和 Sharp [25] 讨论了他们使用该程序的经验。Bhansali [26] 报告了非常令人失望的结果,并声称这些结果是由赤池方法获得。实际上,这些令人失望的结果是由于他对关联统计量的不正确定义造成的,与本文所提出的最小 FPE 程序无关。该程序还扩展到了多变量 AR 模型的拟合 [27]。基于这一识别程序所获得的结果成功应用于水泥窑过程的计算机控制中,其成果由 Otomo 等人报道 [28]。本节讨论的三个程序具有一个共同特性,即对统计量的分析必须扩展到主项 1 / N 1/N 1/N 的阶。
4.以对数似然均值作为拟合优度的度量
众所周知,在正则条件下,MLE 是渐近有效的 [29],这表明似然函数往往是衡量模型参数与真实值微小偏差最为敏感的准则。考虑如下情况: x 1 , x 2 , ⋯ , x N x_{1},x_{2},\cdots,x_{N} x1,x2,⋯,xN 是对一个随机变量进行 N N N 次独立观测,其概率密度函数为 g ( x ) g(x) g(x)。如果给定一个参数族密度函数 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ),其中参数为向量 θ \theta θ,则平均对数似然(或对数似然除以 N N N)定义为
( 1 / N ) ∑ i = 1 N log f ( x i ∣ θ ) (1) (1 / N) \sum_{i=1}^{N} \log f\left(x_{i} \mid \theta\right) \tag{1} (1/N)i=1∑Nlogf(xi∣θ)(1)
(在本文后续部分中,log 表示自然对数。)当 N N N 无限增加时,这个平均值以概率 1 收敛到
S ( g ; f ( ⋅ ∣ θ ) ) = ∫ g ( x ) log f ( x ∣ θ ) d x S(g ; f(\cdot \mid \theta))=\int g(x) \log f(x \mid \theta) \, d x S(g;f(⋅∣θ))=∫g(x)logf(x∣θ)dx
这里假定该积分存在。由 MLE 的有效性可知,对数似然均值 S ( g ; f ( ⋅ ∣ θ ) ) S(g ; f(\cdot \mid \theta)) S(g;f(⋅∣θ)) 必然是对 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 与 g ( x ) g(x) g(x) 之间微小偏离最敏感的准则。其差值
I ( g ; f ( ⋅ ∣ θ ) ) = S ( g ; g ) − S ( g ; f ( ⋅ ∣ θ ) ) I(g ; f(\cdot \mid \theta))=S(g ; g)-S(g ; f(\cdot \mid \theta)) I(g;f(⋅∣θ))=S(g;g)−S(g;f(⋅∣θ))
被称为对 g ( x ) g(x) g(x) 与 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 的 Kullback-Leibler 平均信息,仅当 f ( x ∣ θ ) = g ( x ) f(x \mid \theta)=g(x) f(x∣θ)=g(x)(几乎处处成立)时取零值,否则取正值 [30]。这些观察表明,通过最大化 S ( g ; f ( ⋅ ∣ θ ) ) S(g ; f(\cdot \mid \theta)) S(g;f(⋅∣θ)) 或从熵的概念类比出发最小化 − S ( g ; f ( ⋅ ∣ θ ) ) -S(g ; f(\cdot \mid \theta)) −S(g;f(⋅∣θ)),可以用来定义最佳拟合模型。这里应当提到,早在 1950 年,Bartlett 就采用最后这一量作为信息函数的定义 [31]。 S ( g ; f ( ⋅ ∣ θ ) ) S(g ; f(\cdot \mid \theta)) S(g;f(⋅∣θ)) 最重要的特性之一在于,它的自然估计——即式 (1) 所表示的对数似然均值——可以在不知道 g ( x ) g(x) g(x) 的情况下得到。当只有单一族 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 时,对式 (1) 的极大化关于 θ \theta θ 就得到了 MLE θ ^ \hat{\theta} θ^。
在统计识别问题中,通常会给出几个族的 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ),这些族可能具有不同的形式和/或虽然形式相同但对参数向量 θ \theta θ 的约束不同,需要决定哪一个 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 是最优选择。经典的最大似然原则无法为这类问题提供有用的解决方案。一种解决方案是将前一节中讨论的基本思想与最大似然原则结合起来。
设想如下情况: g ( x ) = f ( x ∣ θ 0 ) g(x)=f\left(x \mid \theta_{0}\right) g(x)=f(x∣θ0)。在这种情形下,我们将 I ( g ; f ( ⋅ ∣ θ ) ) I(g ; f(\cdot \mid \theta)) I(g;f(⋅∣θ)) 与 S ( g ; f ( ⋅ ∣ θ ) ) S(g ; f(\cdot \mid \theta)) S(g;f(⋅∣θ)) 分别记为 I ( θ 0 ; θ ) I\left(\theta_{0} ; \theta\right) I(θ0;θ) 和 S ( θ 0 ; θ ) S\left(\theta_{0} ; \theta\right) S(θ0;θ)。当 θ \theta θ 与 θ 0 \theta_{0} θ0 足够接近时, I ( θ 0 ; θ ) I\left(\theta_{0} ; \theta\right) I(θ0;θ) 可近似为 [30]
I ( θ 0 ; θ 0 + Δ θ ) = ( 1 2 ) ∥ Δ θ ∥ J 2 I\left(\theta_{0} ; \theta_{0}+\Delta \theta\right)=\left(\frac{1}{2}\right)\|\Delta \theta\|_{J^{2}} I(θ0;θ0+Δθ)=(21)∥Δθ∥J2
其中
∥ Δ θ ∥ J 2 = Δ θ ′ J Δ θ \|\Delta \theta\|_{J^{2}}=\Delta \theta^{\prime} J \Delta \theta ∥Δθ∥J2=Δθ′JΔθ
且 J J J 为 Fisher 信息矩阵,其是正定的,定义为
J i j = E { ∂ log f ( X ∣ θ ) ∂ θ i ∂ log f ( X ∣ θ ) ∂ θ j } J_{ij}=E\left\{\frac{\partial \log f(X \mid \theta)}{\partial \theta_{i}} \frac{\partial \log f(X \mid \theta)}{\partial \theta_{j}}\right\} Jij=E{∂θi∂logf(X∣θ)∂θj∂logf(X∣θ)}
其中 J i j J_{ij} Jij 表示矩阵 J J J 的第 ( i , j ) (i, j) (i,j) 个元素, θ i \theta_{i} θi 为 θ \theta θ 的第 i i i 个分量。因此,当 θ 0 \theta_{0} θ0 的 MLE θ ^ \hat{\theta} θ^ 非常接近 θ 0 \theta_{0} θ0 时,由 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 定义的分布与真实分布 f ( x ∣ θ 0 ) f\left(x \mid \theta_{0}\right) f(x∣θ0) 在 S ( g ; f ( ⋅ ∣ θ ) ) S(g ; f(\cdot \mid \theta)) S(g;f(⋅∣θ)) 的变化意义下的偏离就可以用
( 1 2 ) ∥ θ − θ 0 ∥ J 2 \left(\frac{1}{2}\right)\left\|\theta-\theta_{0}\right\|_{J^{2}} (21)∥θ−θ0∥J2
来衡量。现在考虑如下情况:为了最大化似然函数,参数 θ \theta θ 的变化被限制在一个不包含 θ 0 \theta_{0} θ0 的低维子空间 θ \theta θ 内。对于在该子空间内求得的 θ 0 \theta_{0} θ0 的 MLE θ ^ \hat{\theta} θ^,如果使 S ( θ 0 ; θ ) S\left(\theta_{0} ; \theta\right) S(θ0;θ) 取得最大值的 θ \theta θ 足够接近 θ 0 \theta_{0} θ0,则可以证明,在一定正则条件下,当 N N N 足够大时,随机变量
N ∥ θ ^ − θ ∥ J 2 N\|\hat{\theta}-\theta\|^2_J N∥θ^−θ∥J2
的分布近似服从自由度等于该限制参数空间维数的卡方分布(例如参见 [32])。于是有如下关系成立:
E ∞ 2 N I ( θ 0 ; θ ^ ) = N ∥ θ − θ 0 ∥ J 2 + k (2) E_{\infty} 2 N I\left(\theta_{0} ; \hat{\theta}\right)=N\left\|\theta-\theta_{0}\right\|_{J^{2}}+k \tag{2} E∞2NI(θ0;θ^)=N∥θ−θ0∥J2+k(2)
其中 E ∞ E_{\infty} E∞ 表示该近似分布的均值, k k k 是 θ \theta θ 的维数,即为求得 θ ^ \hat{\theta} θ^ 时独立调节参数的个数。关系 (2) 是前一节中所讨论的预测误差期望的推广。当存在多个模型时,自然会选择使 E I ( θ 0 ; θ ^ ) E I\left(\theta_{0} ; \hat{\theta}\right) EI(θ0;θ^) 最小的那个模型。为此,在考虑这些模型的 θ \theta θ 都足够接近 θ 0 \theta_{0} θ0 的情形下,就有必要发展出对 (2) 中 N ∥ θ − θ 0 ∥ J 2 N\|\theta-\theta_{0}\|_{J^{2}} N∥θ−θ0∥J2 的某种估计。关系 (2) 基于如下事实: N ( θ ^ − θ ) \sqrt{N}(\hat{\theta}-\theta) N(θ^−θ) 的渐近分布近似于均值为零、方差矩阵为 J − 1 J^{-1} J−1 的 Gaussian 分布。因此,如果采用
2 ( ∑ i = 1 N log f ( x i ∣ θ 0 ) − ∑ i = 1 N log f ( x i ∣ θ ^ ) ) (3) 2\left(\sum_{i=1}^{N} \log f\left(x_{i} \mid \theta_{0}\right)-\sum_{i=1}^{N} \log f\left(x_{i} \mid \hat{\theta}\right)\right) \tag{3} 2(i=1∑Nlogf(xi∣θ0)−i=1∑Nlogf(xi∣θ^))(3)
作为 N ∥ θ − θ 0 ∥ J 2 N\|\theta-\theta_{0}\|_{J^2} N∥θ−θ0∥J2 的估计,则由于用 θ ^ \hat{\theta} θ^ 替代 θ \theta θ 会引入向下偏差,故需要对 (3) 作修正,即简单地在 (3) 上加上 k k k。对于识别问题来说,仅需要比较各模型下对 E I ( θ 0 ; θ ^ ) E I\left(\theta_{0} ; \hat{\theta}\right) EI(θ0;θ^) 的估计值,因此包含 θ 0 \theta_{0} θ0 的公共项可以省略。
5.信息准则的定义
基于前述观察,引入信息准则 AIC,其定义为
AIC ( θ ^ ) = ( − 2 ) log (最大似然) + 2 k \operatorname{AIC}(\hat{\theta})=(-2) \log \text{(最大似然)}+2k AIC(θ^)=(−2)log(最大似然)+2k
其中,如前所述, k k k 是获得 θ ^ \hat{\theta} θ^ 时独立调节参数的个数。 ( 1 / N ) AIC ( θ ^ ) (1/N)\operatorname{AIC}(\hat{\theta}) (1/N)AIC(θ^) 可看作是 − 2 E S ( θ 0 ; θ ^ ) -2E S\left(\theta_{0} ; \hat{\theta}\right) −2ES(θ0;θ^) 的一个估计。IC 表示信息准则,之所以在前面加上 A,是为了后续可能出现诸如 BIC、DIC 等类似统计量。当存在若干个对应于不同模型的 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 规格时,MAICE 被定义为使得 AIC ( θ ^ ) \operatorname{AIC}(\hat{\theta}) AIC(θ^) 取得最小值的 f ( x ∣ θ ^ ) f(x \mid \hat{\theta}) f(x∣θ^)。而当只有一个无约束族的 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 存在时,则 MAICE 定义为 f ( x ∣ θ ^ ) f(x \mid \hat{\theta}) f(x∣θ^),其中 θ ^ \hat{\theta} θ^ 与经典的 MLE 相同。需要注意的是,当不考虑不同观测集之间结果的比较时,可以在 AIC ( θ ^ ) \operatorname{AIC}(\hat{\theta}) AIC(θ^) 的定义中加入任意的加性常数。当前的 MAICE 定义给出了模型构造中简约原则的数学形式。当两个模型的最大似然值相同时,MAICE 选择参数数较少的模型。
在时间序列分析中,即使在 Gaussian 假定下,似然的精确定义通常也过于复杂而不便实际使用,此时需要作出一定的近似。对于 MAICE 的应用,在如何定义似然函数近似的问题上存在一个微妙问题。这是由于在 AIC 定义中,对数似然必须在数量级 1 1 1 上被一致地定义。对于拟合平稳 Gaussian 过程模型,可以定义模型与真实结构偏离的一种度量为:当观测数 N N N 无限增加时,平均对数似然的极限值。该量与拟合模型所定义的创新过程的对数似然均值相同。故而,将平稳零均值 Gaussian 过程模型拟合到观测序列 y 1 , y 2 , ⋯ , y N y_{1},y_{2},\cdots,y_{N} y1,y2,⋯,yN 的一种自然程序是:首先定义一个原始平稳 Gaussian 模型,其 l l l 阶协方差矩阵 R ( l ) R(l) R(l) 定义为
R ( l ) = ( 1 / N ) ∑ n = 1 N − l y n + l y n ′ , l = 0 , 1 , 2 , ⋯ , N − 1 , = 0 , l = N , N + 1 , ⋯ \begin{aligned} R(l) &= (1/N)\sum_{n=1}^{N-l} y_{n+l} y_{n}^{\prime},\quad l=0,1,2,\cdots,N-1,\\[1mm] &= 0,\quad l=N,N+1,\cdots \end{aligned} R(l)=(1/N)n=1∑N−lyn+lyn′,l=0,1,2,⋯,N−1,=0,l=N,N+1,⋯
然后,通过极大化创新的对数似然均值,或者等价地,如果创新过程协方差矩阵的各元素位于参数集中,则通过最小化创新方差矩阵的对数行列式(这一项乘以 N N N 后用于替代 AIC 定义中的对数似然)来拟合模型。定义 R ( l ) R(l) R(l) 时采用分母 N N N 是为了保持协方差矩阵序列为正定。通过这种以原始模型拟合 Gaussian 模型的程序,在 [33] 中有详细讨论。这自然而然地引出了 Whittle 所发展的 Gaussian 估计的概念。当 y n y_{n} yn 的归一化相关系数的渐近分布与 Gaussian 过程相同时,由这些相关系数构成的统计量的渐近分布也将不依赖于 Gaussian 过程的假设。关于这一点以及为证明当前 AIC 定义而所需相关统计量的渐近行为,在 Whittle 的上述论文中有详细讨论。对于拟合一元 Gaussian AR 模型,用当前 AIC 定义得到的 MAICE 在渐近上与最小 FPE 程序得到的估计相同。
AIC 以及 MAICE 的原始定义最初由作者于 1971 年提出 [3]。一些早期的成功应用成果见于 [3], [35], [36]。
6. 数值例子
在讨论 MAICE 特性之前,本节先演示其实际效用。为了方便读者自行检验结果,文中将 Gaussian AR 模型拟合于 Anderson 所著《时间序列分析》一书中给出的数据上。对 Wold 的三个系列(这些系列是通过二阶 AR 方程人工生成的)拟合了最高达第 50 阶的模型。在两个案例中,MAICE 选出的模型为二阶模型;而在另一个案例中,MAICE 选出的一阶模型是由于生成方程中二阶系数的绝对值相对于其抽样变异性非常小,以致使用 MAICE 得到的一步预测误差方差比采用 MLE 得到的二阶模型更小。对于经典的 Wolfer 太阳黑子数系列( N = 176 N=176 N=176),在拟合了最高达第 35 阶的 AR 模型后,MAICE 选出的是八阶模型。AIC 在二阶时取得局部最小值。在 Beveridge 的小麦价格指数系列( N = 370 N=370 N=370)中,对最高达第 50 阶的 AR 模型的分析中,MAICE 依然选出了八阶模型,而 AIC 在二阶取得了局部最小值,该二阶模型曾被 Sargan 采用。结合 Anderson 对这些序列的讨论,选择八阶模型对这两个序列来说看起来是合理的。
另外,文中报告了两例使用最小 FPE 程序(其渐近上等价于 MAICE)应用的实例。在 Jenkins 和 Watts 的书 [39, 第 5.4.3 节] 中的实例中,所得到的估计与该书作者经过仔细分析后所选的模型一致。在 Whittle 处理的 seiche(涌浪)记录的实例中 [40],最小 FPE 程序清楚表明需要一个非常高阶的 AR 模型。Whittle 在 [41, p. 38] 中讨论了对这一数据集拟合 AR 模型的困难。
该程序还被应用于 Box 和 Jenkins 所著书中给出的系列 E E E 与 F F F。对于系列 E E E( N = 100 N=100 N=100,该系列是 Wolfer 太阳黑子数系列的一部分),作者建议采用二阶或三阶 AR 模型,而在拟合了最高达第 20 阶的 AR 模型后,MAICE 选出了二阶模型。对于系列 F F F( N = 70 N=70 N=70),在拟合了最高达第 10 阶的 AR 模型后,MAICE 选出的也是二阶模型,这与作者的建议是一致的。
为了检验区分 AR 和 MA 模型的能力,作者生成了10个长度为
n
=
1
,
⋯
,
1600
n=1,\cdots,1600
n=1,⋯,1600 的序列,其生成关系为
y
n
=
x
n
+
0.6
x
n
−
1
−
0.1
x
n
−
2
,
y_{n}=x_{n}+0.6\,x_{n-1}-0.1\, x_{n-2},
yn=xn+0.6xn−1−0.1xn−2,
其中
x
n
x_{n}
xn 来源于物理噪声源,并假定为 Gaussian 白噪声。对每个序列的前
N
N
N 个点分别(
N
=
50
,
100
,
200
,
400
,
800
,
1600
N=50,100,200,400,800,1600
N=50,100,200,400,800,1600)拟合了 AR 模型。随着
N
N
N 依次增加,MAICE 选出的 AR 模型的阶数样本平均值分别为
3.1
,
4.1
,
6.5
,
6.8
,
8.2
3.1,4.1,6.5,6.8,8.2
3.1,4.1,6.5,6.8,8.2 和
9.3
9.3
9.3。随后,将一种旨在为 Markov 模型拟合获得初步 MAICE 估计的近似 MAICE 程序(参见 [33])应用于数据。除极少数例外,该近似 MAICE 均为二阶,这对应于 AR-MA 模型中二阶 AR 和一阶 MA 的结构。然后,对
N
=
1600
N=1600
N=1600 的数据分别拟合了二阶和三阶 MA 模型。在拟合到数据的 AR 与 MA 模型中,二阶 MA 模型被选为 MAICE 9 次,而三阶 MA 模型仅选了 1 次。AR 与 MA 模型之间 AIC 最小值的平均差为 7.7,这大致意味着,对于一组
N
=
1600
N=1600
N=1600 的数据,两拟合模型的期望似然比约为 47,更倾向于 MA 模型。
另一测试例子采用了 Gersch 和 Sharp [25] 讨论的实例。生成了 8 个长度为 N = 800 N=800 N=800 的序列,其生成方式为文中所述的 AR-MA 方案。8 个案例中 MAICE 选出的 AR 模型平均阶数为 17.9,这与 Gersch 和 Sharp 报告的值相符。近似 MAICE 程序被用于确定该过程的马尔可夫表示的阶数或维数。对于这 8 个案例,该程序均正确地选择了阶数为 4。对同一数据集,作者拟合了不同阶次的 AR-MA 模型,并计算了相应的 AIC ( p , q ) \operatorname{AIC}(p, q) AIC(p,q) 值,其中 AIC ( p , q ) \operatorname{AIC}(p, q) AIC(p,q) 是指 AR 阶 p p p 与 MA 阶 q q q 模型下的 AIC 值,其定义为
AIC ( p , q ) = N log (MLE of innovation variance) + 2 ( p + q ) \operatorname{AIC}(p, q)=N \log \text{(MLE of innovation variance)}+2(p+q) AIC(p,q)=Nlog(MLE of innovation variance)+2(p+q)
结果为: AIC ( 3 , 2 ) = 192.72 , AIC ( 4 , 3 ) = 66.54 , AIC ( 4 , 4 ) = 67.44 , AIC ( 5 , 3 ) = 67.48 , AIC ( 6 , 3 ) = 67.65 , \operatorname{AIC}(3,2)=192.72,\ \operatorname{AIC}(4,3)=66.54,\ \operatorname{AIC}(4,4)=67.44,\ \operatorname{AIC}(5,3)=67.48,\ \operatorname{AIC}(6,3)=67.65, AIC(3,2)=192.72, AIC(4,3)=66.54, AIC(4,4)=67.44, AIC(5,3)=67.48, AIC(6,3)=67.65, 以及 AIC ( 5 , 4 ) = 69.43 \operatorname{AIC}(5,4)=69.43 AIC(5,4)=69.43。最小值在 p = 4 p=4 p=4 和 q = 3 q=3 q=3 时取得,这对应于真实结构。图 1 展示了通过对该数据集应用各种方法得到的功率谱估计结果。需要说明的是,在此例中,当模型的 p p p 和 q q q 同时大于 4 和 3 时,平均对数似然函数的 Hessian 矩阵在真实参数值处变为奇异。关于这种奇异性问题的详细讨论超出本文讨论范围。图 2 展示了同类方法在一段脑电波记录( N = 1420 N=1420 N=1420)上的应用结果。在此例中,仅拟合了一个 AR-MA 模型,该模型的 AR 阶数为 4,MA 阶数为 3,其 AIC 值为 1145.6,而 MAICE 所选 AR 模型的 AIC 值为 1120.9。这表明 13 阶的 MAICE AR 模型更为合适,这一结论与从图 2 得到的直观印象基本一致。
7.讨论
当 f ( x ∣ θ ) f(x \mid \boldsymbol{\theta}) f(x∣θ) 与 g ( x ) g(x) g(x) 差距非常大时, S ( g ; f ( ⋅ ∣ θ ) ) S(g ; f(\cdot \mid \boldsymbol{\theta})) S(g;f(⋅∣θ)) 仅仅是对 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 与 g ( x ) g(x) g(x) 之间偏差的一种主观度量。因此,对 MAICE 特性的一般讨论只有在至少存在一族 f ( x ∣ θ ) f(x \mid \boldsymbol{\theta}) f(x∣θ) 与 g ( x ) g(x) g(x) 足够接近,从而远小于因用 θ ^ \hat{\theta} θ^ 替代 θ \theta θ 而产生的偏差预期时才有可能。只有当存在几个族满足这一条件时,才有必要对 MAICE 的统计特性进行详细分析。作为对 − 2 N E S ( g ; f ( ⋅ ∣ θ ^ ) ) -2NE S(g ; f(\cdot \mid \hat{\theta})) −2NES(g;f(⋅∣θ^)) 的一个单一估计, − 2 -2 −2 乘以对数最大似然即可满足要求,但为达到“估计两者差异”的目的,在 AIC 的定义中引入加项 + 2 k +2k +2k 则尤为关键。Bhansali [26] 所报道的令人失望的结果正是由于其错误地使用了等价于在 AIC 中使用 + k +k +k 而非 + 2 k +2k +2k 的统计量。
当模型通过对 f ( x ∣ θ ) f(x \mid \theta) f(x∣θ) 的参数 θ \theta θ施加逐步限制而得到时,MAICE 程序实际上表现为对传统对数似然比拟合优度检验的反复应用,其显著性水平均由加项 + 2 k +2k +2k 自动调整。当存在不同族在逼近真实似然方面同样良好时,局部上情况可以近似为同一族的不同参数化。对这些情况,两模型间 AIC 差值的显著性可通过将它与一个自由度等于两模型参数 k k k 之差的卡方变量的变异性比较来评定。当两模型在 Cox [42], [43] 定义意义上构成彼此独立的不同族时,Cox 所发展并由 Walker [44] 扩展到时间序列情况的程序可能有助于对 AIC 差值进行详细评估。
必须明确认识到,除非将假设检验程序定义为具有所需显著性水平的决策程序,否则 MAICE 无法与之比较。对于具有不同比例参数数目的模型,使用固定的显著性水平进行比较是不正确的,因为这并未考虑到参数数目增加时估计变异性的增加。正如 Kennedy 和 Bancroft [45] 的工作所显示,基于一系列显著性检验所构建的模型构建理论尚未充分发展到能够提供实际有用程序的程度。
尽管作者并未证明 MAICE 的最优性,目前它是唯一适用于能够正确定义似然的所有情形的程序,而且在实际中能够产生非常合理的结果,而几乎无需主观判断。数值实验的成功结果表明,MAICE 在建模、预测、信号检测、模式识别和自适应等领域具有近乎无限的应用潜力。对于 AIC 定义和使用方法的进一步改进,以及在各具体应用中 MAICE 同其它程序之间的数值比较,将成为今后研究的主题。
8.结论
作为统计模型构造或识别方法,假设检验程序的实际效用应被认为相当有限。为了发展有用的识别程序,有必要采取一种更直接控制由所选模型所引起的误差或损失的方法。鉴于经典最大似然程序的成功,平均对数似然看似是衡量统计模型拟合优度的一个自然选择。基于 AIC 的 MAICE 程序不仅提供了一种多用途的统计模型识别方法,而且也为模型构造中简约原则提供了数学表述。由于基于 MAICE 的程序的实现几乎不需要主观判断,数值应用上的成功结果表明,许多用于预测、信号检测、模式识别和自适应等的统计识别程序将有可能借助 MAICE 而变得切实可用。
9.致谢
作者感谢 Stanford University 的 Kailath 教授鼓励其撰写本论文。同时,亦感谢长崎大学的 Sato 教授提供了文中第 V 节所使用的脑电波数据。