时间序列分析与动态数据建模

构造和式
由于 在单位圆内没有零点，上式被积函数除了r=0时有极点在原点外，r≥1时在单位圆内是解析的，根据柯西留数定理可得
详细写，就是
若令
(5-2-7)
则有
(5-2-8)
利用式（5-2-7）可将式（5-2-6）写成
(5-2-9)
其中
将式（5-2-9）代入式（5-2-2），并考虑到 的偶函数性质，得
(5-2-10)
由结果可以看出，在已知头M+1个相关函数值时，以它作为约束条件推出的极大熵谱是AR模型的功率谱。
一般情况下，我们直接观察到的是过程 的数据，并不知道相关函数的准确值。因此通常根据 为最小来求 ，这将得到一组n个方程，其形式和式（5-2-8）一样，只是用估计的自相关代替R(r)，所以当采用这种估计值时，极大熵法和最小二乘法估计的结果是相同的。
阶数从1起 和 的递推计算将利用自相关函数阵的对称托普里茨（Toeplitz）性质，以二阶向三阶递推为例，由式（5-2-8）写出
(5-3-2)
如果将该矩阵所对应的议程组及变量的顺序反过来，则有
(5-3-3)
把式（5-3-2）的相关阵放大到4×4；可得
(5-3-4)
其中
(5-3-5)
对式（5-3-3）也作类似扩大，然后将两个4×4阵按下面方式组合
(5-1-15)
。式（5-1-15）是随机变量 通过线性变换Ax成为随机变量 ，这里
根据式(5-1-12)得
除以（t+1）并令 则
现在只要证明 等于式（5-1-13）中的第二项积分就够了。由于 ， 情况下有
这里的线积分是沿单位园进行的，因
故式（5-1-13）中的第二项积分等于 ，所以需要证明
由于G(B)在单位园内是解析的，所以上式中的积分路线可以任意小，当 。故上式右边等于 ，证明完成。
5.3 极大熵谱估计的伯格算法
上面已经指出，由于不确切知道头M+1个自相关的真值一般只能用它的估计值代替，在第三章中提到的两种估计算式为
和
是非负定的，方差较小，但估计偏度随r增加而增大。R是渐近无偏的，但不能保证非负定性质。因此以上两种算法各有不足之处。伯格采用的算法是从一阶模型开始逐步增加阶数的递推算法，每步递推都能保证相应的自相关序列是非负定的，而且得到的模型也是平稳的，不仅如此，从后面介绍还可看出，由于采用了双向（正、反）预报误差平方和为最小，提高了数据的利用率，充分挖取数据内含的信息，因此特别有利于短数据的分析和建模。
（3）若 是正态过程，其功率谱 满足
(5-1-15)
则有
(5-1-16)
注：
这一结论是不难看出的，因为非白正态过程的功率谱密度 可以看作是方差为1的白噪声通过频率响应模的平方等 的线性系统所产生的过程的谱，因此利用式（5-1-13）和（5-1-9）就可导出式（5-1-16）。
式（5-1-16）给出了过程的功率谱密度和它的熵率之间的关系式，由于右边第一项是常数，比较 的大小等价于比较第二项积分的大小，因此称 的谱熵，并以它作为推导极大熵谱估计的出发点。
最小残差方差为0.0304。
图5-3为极大熵谱曲线，其峰值出现在 处。
这里需要指出的是，用自回归模型拟合只适用于有谱密度的序列，对正弦信号而言其谱密度在给定频率处为无穷。这个例子只是说明根据短的序列样本以极大熵谱估计谱线的位置，即正弦振荡的频率。
极大熵谱估计的伯格算法程序见附录二MEBURG。
5.4极大熵谱估计的LS—LUD算法
下面的例子是一组由
产生的20个数据 其中 是白噪声，它的标准差约为正弦振幅的5%，一个 的纪录为0.1410，1.0509，1.7826。2.6804，3.0536, 2.9605, 2.7524, 2.1767, 1.6413, 1.371, 0.1217, -0.9359, -1.8501, -2.5495, -2.5454, -2.9358, -3.0448, -2.2961, -1.7726, -0.9091（见图5-2），利用MEBURG程序计算可得最佳阶数
其中A是N×N非奇异矩阵，X的联合概率密度为 ，则由于
(5-1-11)
可得
(5-1-12)
（2）若 是一个稳定的因果系统的输入，该系统的传递函数为G(B)（这在单位园内无极点），系统单位脉冲响应为 。设 在 时开始输入，因而系统输出 是平稳的，以
(5-1-13)
其中
(5-1-14)
证明：若 在t=0时开始输入，则系统输出为
例如已知过程 的方差为 ，即
(5-1-17)
要导出能使 为最大的功率谱 。这个问题可以通过求泛函极值来解决。以 表拉格伦日乘子作Hale Waihona Puke Baidu函
其变为
达到极值的条件为 ,故应有 代回式（5-1-17）可得 ，即 必须为常数 ，因此只有当过程为白噪声时才能使熵率达到最大，这里约束条件是方差为 。
5.2 极大熵准则的谱估计
式（5-1-6）就是长度为N的正态时间序列的熵。若有正态白噪声 （方差为 ），则
可求得其熵为
(5-1-7)
由于H随N增长而发散，定义熵率h为
(5-1-8)
故白噪声过程的熵率为
(5-1-9)
3．功率谱和熵率的关系
下面给出功率谱和熵率间的一些重要性质和关系。
（1）如果随机向量 是随机向量 ，则由于
(5-1-10)
由 得
(5-2-2)
这里应有
(5-2-3)
以保证 是实的。将式（5-2-2）代入式（5-2-1）得
(5-2-4)
用后移算子 代替变量 ，上式所写成
(5-2-5)
该积分沿B平面单位园进行，基于式（5-2-3）有
(5-2-6)
其中：
不难看出
故
多项式 的全部零点均在B平面单位圆外，而 的全部零点均在单位圆内，两部分零点是互为倒数分布的。将式（5-2-6）代入式（5-2-5）得
第五章目录
第五章 极大熵谱估计
1967年伯格(J．P．Burg)刚一发表：极大熵谱分析”的方法就在工程和科技界产生很大影响，成为相当流行的功率谱密度估计方法。伯格在谱估计准则的提出和具体算法上有所创新，由此演变出来的算法有很多种，被统称为“现代谱分析”。
5.1谱熵和极大熵准则
1．问题的提出
从19世纪未舒斯特(Schuster)在利用富氏级数分析信号隐含的周期特性时提出了“周期图”，到1985年由伯来克曼和杜奇提出了谱估计的“间接法”和1965年FFT算法提出后流行的“直接法”，它们本质上都是把原序列经过开窗截取处理来获得对序列谱密度的估计。不论对数据加窗还是对自相关函数加窗，其目的都在于使谱估计的方差减小，然而加窗不可避免地产生频域“泄漏”，使功率谱失真，尽管在窗函数形式的选择和处理方法上做了很多分析研究，使得以周期图为基础的方法达到相当成熟和实用的程度，但是任何抑制旁瓣的方法都是以损失谱分辨力为代价的，这个难题在数据量少的情况下更为突出。
2．高斯过程的熵和熵率
假定我们研究的随机试验a只有有限个不相容的结果 ，它们相应的概率为 ，且满足 ，简单描述如下：
申农找到并证明了可以用 这个量来度量 的不肯定性的程度：
或简写成：
称为试验 的熵
当随机变量的可能取值是连续的，则H定义式中的和式用积分代替
(5-1-1)
其中 为随机变量，对数可以取10或取e为底，在比较熵的大小时并没有影响，下面为计算方便均以自然对数ln来定义，如x为正态随机变量， ，则有
(5-1-2)
进一步，如果讨论的是时间序列的实现 则这一过程的熵用下面N维积分表示：
(5-1-3)
其中 是联合概率密度函数 ，若时间序列是高斯的，则
(5-1-4)
其中 为自协方差阵
(5-1-5)
它的i行j列元素为 的均值， 表均值向量。
将式（5-1-4）代入式（5-1-3）求过程x的熵
(5-1-6)
数学家申农最早提出“熵”的概念，在统计学中用它作为各种随机试验的不肯定性程度的度量。在热力学和信息论中，“熵”都有其具体的物理背景和应用。后面介绍将会看到，满足熵极大的谱估计是自回归模型的谱。1971年凡登包士（Van Den Bos）证明，一维极大熵谱估计和自回归谱的最小二乘估计是等效的。尽管如此，伯格关于熵谱估计的概念和他对自回归参数的递推算法却独树一帜，随后还有人提出了各种改进算法，但要注意把极大熵概念本身同等法区别开来。
对它作极小化求
然而伯格指出，这样求得的模型参数并不总能满足平稳性条件，但他注意到利文森（Levinson）提出的递推算法可以做到这一点，这种算法由AR(1)推出AR（2）是按以下关系：
伯格决定不由 ，作极小化，而是把 看作是 的函数，然后对 极小化 ，为此写出
其中
表示一阶模型的正向预报误差， 表示相应的反向预报误差。由 可得
（1）递推所得的参数满足平稳性条件，即
的根全部在B平面单位圆以外，或者等效地说
中任一 都满足 。
证：令 并将 表为 ，则 B应当是一个半径小于1的圆，或者说 是圆心在(1,0)但不包围原点的圆（见图5-1）。这必等价于当f由-1/2时变到
1/2时 的幅角增量为零。当全部 的模均小于1时， 总的幅角增量也是零，或 曲线不包围原点。
先看一阶模型
这里上标是用来标明它属于一阶模型，以下都用它作为递推次数的记号，对平方和
作极小化来选择 时可能会出现 ，从而使模型不平稳，例如序列值为1，2时
由
由于在均方意义下最优的模型参数只取决于序我的自相关而非序列值本身，而序列按反向的时间顺序排列并不改变自相关函数（见2.1节），伯格提出以
作极小化来求 ，上式右边第二个括弧内的 可以看作是由 “预报” 的误差，不难看出，由此得到的 满足 。因此这种将正反双方预报误差的平方和作极小化的办法在一阶的情况下是可行的如果自然地推广到二阶，则有
根据伯格所提出的概念，功纺谱密度估计的准则应当是：
设 表示估计的谱，则它在满足约束条件
(5-2-1)
的同时，应使谱熵 达到极大，其中 是 的正、实、偶函数，这样对应的R(r)自然也是r的偶函数。
下面论证满足以上要求的 所应具有的形式。
设已知自相关函数R(r)在 内的2M＋1个值，以 表拉格伦日乘子作泛函
问题的实质是：在周期图估计中，我们对数据或是它的相关函数所做的加窗处理，等于是假定在窗口外数据(或自相关)为零，而窗口内的部分则加上某种形式的修正。这些人为措施使来自观察的信息受到了一定程度的歪曲。
伯格提出的新概念是；和估计的功率谱相对应的自相关和由观察数据算得的自相关一致，同时对已有的区段之外的自相关值采用外推的办法求取，而不是一概假定为零，外推的原则是使相应的序列在未知点上取值的可能性具有最大的不确定性，亦即不对结果人为地强添任何增加的信息。
证：由于 ，根据式（5-3-9）必有 。再以 表示由 ， …， 构成的相关函数阵，则式（5-2-8）可写成
由克莱姆法则知
由于 按归纳法可知每次递推所得的 ，因此上述递推算法得出的 序列构成正定列。
关于由极大熵谱获得的模型阶数问题，由式（5-2-10）可见，其阶数M是已给自相关估值的最大迟后，当数据个数为N时最大可能的迟后值为N-1，这可能并非是过程的真正阶数。而另一方面，如果序列本身是无限阶的AR模型（如ARMA模型的等效），需要很高的阶数才能逼远真正的过程，这时已给相关的最大迟后所定出的阶数又可能太小。当然，M愈大，用 估计 的精度也愈低，所以取很大的阶数未必就好。鉴于极大熵是AR谱，我们可以利用诸如FPE、AIC、BIC等定阶准则进行检验和判定。
因 恒大于零，可以证明 。后面将会看到这对模型的平稳性和自相关函数的非负定性都是很关键的。
求 后，AR(2)的参数是
从AR(2)向AR(3)递推的方式和前面类似，令
于是
其中 分别为二阶模型的正掺向预报误差
同样，由 求得
且有
这样的梯推继续下去，到M阶时的一般算式有
(5-3-1)
现在来看 和自相关函数值的递推式，令
(5-3-6)
根据式（5-3-1）中关于自回归参数的递推关系可见上式左边有关参数的矩阵乃是三阶模型的参数 ，因此上式右边等于 ，于是有
(5-3-7)
解之得
利用式（5-3-5）和式（5-3-7）得
推广到一般，可综合如下递推算式：
(5-3-8)
(5-3-9)
(5-3-10)
以上三式和（5-3-1）诸式组成了一套极大熵谱估计的伯格递推算法（程序见附录二MEBURG）。
而今
(5-3-11)
这里模
故当
是不包围原点的，即式（5-3-11）右边[·]部分的零点都在B平面单位圆外，如果前一步递推得到的 已满足平稳条件，则 也将满足平稳条件。由于从一阶开始递推时已有 ，且以后每步递推均有 ，因此每步递推所得的参数必然均能满足平稳条件。
（2）递推所得的自相关序列满足非负定条件。