時間序列分析與動態(tài)數(shù)據(jù)建模

1、第五章目錄第五章 極大熵譜估計(jì)15.1 譜熵和極大熵準(zhǔn)則11問題的提出12高斯過程的熵和熵率13功率譜和熵率的關(guān)系35.2 極大熵準(zhǔn)則的譜估計(jì)65.3 極大熵譜估計(jì)的伯格算法95.4 極大熵譜估計(jì)的LSLUD算法16第五章 極大熵譜估計(jì)1967年伯格(JPBurg)剛一發(fā)表：極大熵譜分析”的方法就在工程和科技界產(chǎn)生很大影響，成為相當(dāng)流行的功率譜密度估計(jì)方法。伯格在譜估計(jì)準(zhǔn)則的提出和具體算法上有所創(chuàng)新，由此演變出來的算法有很多種，被統(tǒng)稱為“現(xiàn)代譜分析”。5.1 譜熵和極大熵準(zhǔn)則1問題的提出從19世紀(jì)未舒斯特(Schuster)在利用富氏級數(shù)分析信號隱含的周期特性時提出了“周期圖”，到1985年由

2、伯來克曼和杜奇提出了譜估計(jì)的“間接法”和1965年FFT算法提出后流行的“直接法”，它們本質(zhì)上都是把原序列經(jīng)過開窗截取處理來獲得對序列譜密度的估計(jì)。不論對數(shù)據(jù)加窗還是對自相關(guān)函數(shù)加窗，其目的都在于使譜估計(jì)的方差減小，然而加窗不可避免地產(chǎn)生頻域“泄漏”，使功率譜失真，盡管在窗函數(shù)形式的選擇和處理方法上做了很多分析研究，使得以周期圖為基礎(chǔ)的方法達(dá)到相當(dāng)成熟和實(shí)用的程度，但是任何抑制旁瓣的方法都是以損失譜分辨力為代價的，這個難題在數(shù)據(jù)量少的情況下更為突出。問題的實(shí)質(zhì)是：在周期圖估計(jì)中，我們對數(shù)據(jù)或是它的相關(guān)函數(shù)所做的加窗處理，等于是假定在窗口外數(shù)據(jù)(或自相關(guān))為零，而窗口內(nèi)的部分則加上某種形式的修正

3、。這些人為措施使來自觀察的信息受到了一定程度的歪曲。伯格提出的新概念是；和估計(jì)的功率譜相對應(yīng)的自相關(guān)和由觀察數(shù)據(jù)算得的自相關(guān)一致，同時對已有的區(qū)段之外的自相關(guān)值采用外推的辦法求取，而不是一概假定為零，外推的原則是使相應(yīng)的序列在未知點(diǎn)上取值的可能性具有最大的不確定性，亦即不對結(jié)果人為地強(qiáng)添任何增加的信息。數(shù)學(xué)家申農(nóng)最早提出“熵”的概念，在統(tǒng)計(jì)學(xué)中用它作為各種隨機(jī)試驗(yàn)的不肯定性程度的度量。在熱力學(xué)和信息論中，“熵”都有其具體的物理背景和應(yīng)用。后面介紹將會看到，滿足熵極大的譜估計(jì)是自回歸模型的譜。1971年凡登包士（Van Den Bos）證明，一維極大熵譜估計(jì)和自回歸譜的最小二乘估計(jì)是等效的。盡管

4、如此，伯格關(guān)于熵譜估計(jì)的概念和他對自回歸參數(shù)的遞推算法卻獨(dú)樹一幟，隨后還有人提出了各種改進(jìn)算法，但要注意把極大熵概念本身同等法區(qū)別開來。2高斯過程的熵和熵率 假定我們研究的隨機(jī)試驗(yàn)a只有有限個不相容的結(jié)果，它們相應(yīng)的概率為，且滿足，簡單描述如下：申農(nóng)找到并證明了可以用這個量來度量的不肯定性的程度：或簡寫成：稱為試驗(yàn)的熵當(dāng)隨機(jī)變量的可能取值是連續(xù)的，則H定義式中的和式用積分代替 (5-1-1)其中為隨機(jī)變量，對數(shù)可以取10或取e為底，在比較熵的大小時并沒有影響，下面為計(jì)算方便均以自然對數(shù)ln來定義，如x為正態(tài)隨機(jī)變量，則有 (5-1-2)進(jìn)一步，如果討論的是時間序列的實(shí)現(xiàn)則這一過程的熵用下面N維

5、積分表示： (5-1-3)其中是聯(lián)合概率密度函數(shù)，若時間序列是高斯的，則 (5-1-4)其中為自協(xié)方差陣 (5-1-5)它的i行j列元素為的均值，表均值向量。將式（5-1-4）代入式（5-1-3）求過程x的熵 (5-1-6)式（5-1-6）就是長度為N的正態(tài)時間序列的熵。若有正態(tài)白噪聲（方差為），則可求得其熵為 (5-1-7)由于H隨N增長而發(fā)散，定義熵率h為 (5-1-8)故白噪聲過程的熵率為 (5-1-9)3功率譜和熵率的關(guān)系 下面給出功率譜和熵率間的一些重要性質(zhì)和關(guān)系。（1）如果隨機(jī)向量是隨機(jī)向量，則由于 (5-1-10)其中A是N×N非奇異矩陣，X的聯(lián)合概率密度為，則由于 (

6、5-1-11)可得 (5-1-12)（2）若是一個穩(wěn)定的因果系統(tǒng)的輸入，該系統(tǒng)的傳遞函數(shù)為G(B)（這在單位園內(nèi)無極點(diǎn)），系統(tǒng)單位脈沖響應(yīng)為。設(shè)在時開始輸入，因而系統(tǒng)輸出是平穩(wěn)的，以 (5-1-13)其中 (5-1-14)證明：若在t=0時開始輸入，則系統(tǒng)輸出為 (5-1-15)。式（5-1-15）是隨機(jī)變量通過線性變換Ax成為隨機(jī)變量，這里根據(jù)式(5-1-12)得除以（t+1）并令則現(xiàn)在只要證明等于式（5-1-13）中的第二項(xiàng)積分就夠了。由于，情況下有這里的線積分是沿單位園進(jìn)行的，因故式（5-1-13）中的第二項(xiàng)積分等于，所以需要證明由于G(B)在單位園內(nèi)是解析的，所以上式中的積分路線可以任

7、意小，當(dāng)。故上式右邊等于，證明完成。（3）若是正態(tài)過程，其功率譜滿足 (5-1-15)則有 (5-1-16)注：這一結(jié)論是不難看出的，因?yàn)榉前渍龖B(tài)過程的功率譜密度可以看作是方差為1的白噪聲通過頻率響應(yīng)模的平方等的線性系統(tǒng)所產(chǎn)生的過程的譜，因此利用式（5-1-13）和（5-1-9）就可導(dǎo)出式（5-1-16）。式（5-1-16）給出了過程的功率譜密度和它的熵率之間的關(guān)系式，由于右邊第一項(xiàng)是常數(shù)，比較的大小等價于比較第二項(xiàng)積分的大小，因此稱的譜熵，并以它作為推導(dǎo)極大熵譜估計(jì)的出發(fā)點(diǎn)。例如已知過程的方差為，即 (5-1-17)要導(dǎo)出能使為最大的功率譜。這個問題可以通過求泛函極值來解決。以表拉格倫日乘子

8、作泛函其變?yōu)檫_(dá)到極值的條件為,故應(yīng)有代回式（5-1-17）可得，即必須為常數(shù)，因此只有當(dāng)過程為白噪聲時才能使熵率達(dá)到最大，這里約束條件是方差為。5.2 極大熵準(zhǔn)則的譜估計(jì)根據(jù)伯格所提出的概念，功紡譜密度估計(jì)的準(zhǔn)則應(yīng)當(dāng)是：設(shè)表示估計(jì)的譜，則它在滿足約束條件 (5-2-1)的同時，應(yīng)使譜熵達(dá)到極大，其中是的正、實(shí)、偶函數(shù)，這樣對應(yīng)的R(r)自然也是r的偶函數(shù)。下面論證滿足以上要求的所應(yīng)具有的形式。設(shè)已知自相關(guān)函數(shù)R(r)在內(nèi)的2M1個值，以表拉格倫日乘子作泛函由得 (5-2-2)這里應(yīng)有 (5-2-3)以保證是實(shí)的。將式（5-2-2）代入式（5-2-1）得 (5-2-4)用后移算子代替變量，上式所

9、寫成(5-2-5)該積分沿B平面單位園進(jìn)行，基于式（5-2-3）有 (5-2-6)其中：不難看出故多項(xiàng)式的全部零點(diǎn)均在B平面單位圓外，而的全部零點(diǎn)均在單位圓內(nèi)，兩部分零點(diǎn)是互為倒數(shù)分布的。將式（5-2-6）代入式（5-2-5）得構(gòu)造和式由于在單位圓內(nèi)沒有零點(diǎn)，上式被積函數(shù)除了r=0時有極點(diǎn)在原點(diǎn)外，r1時在單位圓內(nèi)是解析的，根據(jù)柯西留數(shù)定理可得詳細(xì)寫，就是若令 (5-2-7)則有 (5-2-8)利用式（5-2-7）可將式（5-2-6）寫成 (5-2-9)其中將式（5-2-9）代入式（5-2-2），并考慮到的偶函數(shù)性質(zhì)，得 (5-2-10)由結(jié)果可以看出，在已知頭M+1個相關(guān)函數(shù)值時，以它作為約

10、束條件推出的極大熵譜是AR模型的功率譜。一般情況下，我們直接觀察到的是過程的數(shù)據(jù)，并不知道相關(guān)函數(shù)的準(zhǔn)確值。因此通常根據(jù)為最小來求，這將得到一組n個方程，其形式和式（5-2-8）一樣，只是用估計(jì)的自相關(guān)代替R(r)，所以當(dāng)采用這種估計(jì)值時，極大熵法和最小二乘法估計(jì)的結(jié)果是相同的。5.3 極大熵譜估計(jì)的伯格算法上面已經(jīng)指出，由于不確切知道頭M+1個自相關(guān)的真值一般只能用它的估計(jì)值代替，在第三章中提到的兩種估計(jì)算式為和是非負(fù)定的，方差較小，但估計(jì)偏度隨r增加而增大。R是漸近無偏的，但不能保證非負(fù)定性質(zhì)。因此以上兩種算法各有不足之處。伯格采用的算法是從一階模型開始逐步增加階數(shù)的遞推算法，每步遞推都能

11、保證相應(yīng)的自相關(guān)序列是非負(fù)定的，而且得到的模型也是平穩(wěn)的，不僅如此，從后面介紹還可看出，由于采用了雙向（正、反）預(yù)報(bào)誤差平方和為最小，提高了數(shù)據(jù)的利用率，充分挖取數(shù)據(jù)內(nèi)含的信息，因此特別有利于短數(shù)據(jù)的分析和建模。先看一階模型這里上標(biāo)是用來標(biāo)明它屬于一階模型，以下都用它作為遞推次數(shù)的記號，對平方和作極小化來選擇時可能會出現(xiàn)，從而使模型不平穩(wěn)，例如序列值為1，2時由由于在均方意義下最優(yōu)的模型參數(shù)只取決于序我的自相關(guān)而非序列值本身，而序列按反向的時間順序排列并不改變自相關(guān)函數(shù)（見2.1節(jié)），伯格提出以作極小化來求，上式右邊第二個括弧內(nèi)的可以看作是由“預(yù)報(bào)” 的誤差，不難看出，由此得到的滿足。因此這種

12、將正反雙方預(yù)報(bào)誤差的平方和作極小化的辦法在一階的情況下是可行的如果自然地推廣到二階，則有對它作極小化求然而伯格指出，這樣求得的模型參數(shù)并不總能滿足平穩(wěn)性條件，但他注意到利文森（Levinson）提出的遞推算法可以做到這一點(diǎn)，這種算法由AR(1)推出AR（2）是按以下關(guān)系：伯格決定不由，作極小化，而是把看作是的函數(shù)，然后對極小化，為此寫出其中表示一階模型的正向預(yù)報(bào)誤差，表示相應(yīng)的反向預(yù)報(bào)誤差。由可得因恒大于零，可以證明。后面將會看到這對模型的平穩(wěn)性和自相關(guān)函數(shù)的非負(fù)定性都是很關(guān)鍵的。求后，AR(2)的參數(shù)是從AR(2)向AR(3)遞推的方式和前面類似，令于是其中分別為二階模型的正摻向預(yù)報(bào)誤差同樣

13、，由求得且有這樣的梯推繼續(xù)下去，到M階時的一般算式有 (5-3-1)現(xiàn)在來看和自相關(guān)函數(shù)值的遞推式，令階數(shù)從1起和的遞推計(jì)算將利用自相關(guān)函數(shù)陣的對稱托普里茨（Toeplitz）性質(zhì)，以二階向三階遞推為例，由式（5-2-8）寫出 (5-3-2)如果將該矩陣所對應(yīng)的議程組及變量的順序反過來，則有 (5-3-3)把式（5-3-2）的相關(guān)陣放大到4×4；可得 (5-3-4)其中 (5-3-5)對式（5-3-3）也作類似擴(kuò)大，然后將兩個4×4陣按下面方式組合 (5-3-6)根據(jù)式（5-3-1）中關(guān)于自回歸參數(shù)的遞推關(guān)系可見上式左邊有關(guān)參數(shù)的矩陣乃是三階模型的參數(shù)，因此上式右邊等于，于

14、是有 (5-3-7)解之得利用式（5-3-5）和式（5-3-7）得推廣到一般，可綜合如下遞推算式： (5-3-8) (5-3-9) (5-3-10)以上三式和（5-3-1）諸式組成了一套極大熵譜估計(jì)的伯格遞推算法（程序見附錄二MEBURG）。（1） 遞推所得的參數(shù)滿足平穩(wěn)性條件，即的根全部在B平面單位圓以外，或者等效地說中任一都滿足。證：令并將表為，則B應(yīng)當(dāng)是一個半徑小于1的圓，或者說是圓心在(1,0)但不包圍原點(diǎn)的圓（見圖5-1）。這必等價于當(dāng)f由-1/2時變到圖5-1 的根滿足平穩(wěn)條件時的曲線1/2時的幅角增量為零。當(dāng)全部的模均小于1時，總的幅角增量也是零，或曲線不包圍原點(diǎn)。而今 (5-3

15、-11)這里模故當(dāng)是不包圍原點(diǎn)的，即式（5-3-11）右邊·部分的零點(diǎn)都在B平面單位圓外，如果前一步遞推得到的已滿足平穩(wěn)條件，則也將滿足平穩(wěn)條件。由于從一階開始遞推時已有，且以后每步遞推均有，因此每步遞推所得的參數(shù)必然均能滿足平穩(wěn)條件。（2）遞推所得的自相關(guān)序列滿足非負(fù)定條件。證：由于，根據(jù)式（5-3-9）必有。再以表示由，構(gòu)成的相關(guān)函數(shù)陣，則式（5-2-8）可寫成由克萊姆法則知由于按歸納法可知每次遞推所得的，因此上述遞推算法得出的序列構(gòu)成正定列。關(guān)于由極大熵譜獲得的模型階數(shù)問題，由式（5-2-10）可見，其階數(shù)M是已給自相關(guān)估值的最大遲后，當(dāng)數(shù)據(jù)個數(shù)為N時最大可能的遲后值為N-1，

16、這可能并非是過程的真正階數(shù)。而另一方面，如果序列本身是無限階的AR模型（如ARMA模型的等效），需要很高的階數(shù)才能逼遠(yuǎn)真正的過程，這時已給相關(guān)的最大遲后所定出的階數(shù)又可能太小。當(dāng)然，M愈大，用估計(jì)的精度也愈低，所以取很大的階數(shù)未必就好。鑒于極大熵是AR譜，我們可以利用諸如FPE、AIC、BIC等定階準(zhǔn)則進(jìn)行檢驗(yàn)和判定。下面的例子是一組由產(chǎn)生的20個數(shù)據(jù)其中是白噪聲，它的標(biāo)準(zhǔn)差約為正弦振幅的5%，一個的紀(jì)錄為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.

17、8501, -2.5495, -2.5454, -2.9358, -3.0448, -2.2961, -1.7726, -0.9091（見圖5-2），利用MEBURG程序計(jì)算可得最佳階數(shù)最小殘差方差為0.0304。圖5-3為極大熵譜曲線，其峰值出現(xiàn)在處。圖5-2 20點(diǎn)數(shù)據(jù)圖圖5-3 20點(diǎn)數(shù)據(jù)的極大熵譜AR（6）這里需要指出的是，用自回歸模型擬合只適用于有譜密度的序列，對正弦信號而言其譜密度在給定頻率處為無窮。這個例子只是說明根據(jù)短的序列樣本以極大熵譜估計(jì)譜線的位置，即正弦振蕩的頻率。 極大熵譜估計(jì)的伯格算法程序見附錄二MEBURG。5.4 極大熵譜估計(jì)的LSLUD算法 伯格在極大熵基礎(chǔ)上提

18、出的算法，由于采用了較合理的外推和正反向誤差平方和的極小化，比傳統(tǒng)的方法提高了分辨力。或者說在同樣的分辨力下只需要較少的數(shù)據(jù)量。在伯格方法引起人們廣泛重視和應(yīng)用的同時，也發(fā)現(xiàn)其不足之處，這就是“譜峰偏移”和“譜線分裂”，前者是指峰值頻率估值和真值之間的偏離度，后者是指本來只有一個譜峰，但在估計(jì)譜中卻出現(xiàn)兩個或多個相距很近的譜峰。這類現(xiàn)象在耶爾一瓦克爾的AR譜估計(jì)中也是存在的(凱伊(Kay)和馬波(Marple)已曾指出過)，而伯格算法依然存在這兩個問題。泛杰爾(Fougere)首先指出；當(dāng)數(shù)據(jù)中信噪比高以及所取階數(shù)較高的情況下，伯格算法容易產(chǎn)生譜線分裂，在用周期為了的正弦信號疊加白噪聲作樣本進(jìn)

19、行分析中發(fā)現(xiàn)，當(dāng)數(shù)據(jù)長度為T4的奇數(shù)倍，以及正弦的初始相位為的奇數(shù)倍時也容易出現(xiàn)分裂。譜峰偏移也和初始相位有關(guān)。而數(shù)據(jù)量加大，譜線分裂和譜峰偏移量都會減弱。為了解決上述問題，人們已經(jīng)并且還在做出努力，通常在研究中都采用正反向預(yù)報(bào)誤差，而實(shí)踐結(jié)果表明，AR模型的最小二乘(LS)解的頻率偏移小，在短數(shù)據(jù)下，不受托普里茨(Toeplitz)矩陣形式限制的LS方法可以得到更好的結(jié)果。不過LS方法計(jì)算量大模，型結(jié)果可能非平穩(wěn)等問題則需要妥善解決。這方面近幾年來已取得一些肯定的結(jié)果，本節(jié)介紹的LS-LUD算法就是在最小二乘方法基礎(chǔ)上采用上下三角陣分解(LUD)的算法求解AR譜。我們知道。若給定時間序列為，

20、當(dāng)用n個既往觀察數(shù)據(jù)作一步正向預(yù)報(bào)時。其預(yù)報(bào)誤差為 (5-4-1)其中 (5-4-2)根據(jù)LS準(zhǔn)則，。類似地，當(dāng)用n個后來觀察數(shù)據(jù)作反向預(yù)報(bào)時。其預(yù)報(bào)誤差記作 (5-4-3)其中 (5-4-4)而可由的準(zhǔn)則確定利用矩陣方程表示時。式（5-4-1）可寫成 (5-4-5)其中式（5-4-5）是含有n個未知參數(shù)的N-n個方程組。殘差平方和記作根據(jù)可得n個方程式正規(guī)方程： (5-4-6)或 (5-4-7)其中 類似地，由反向預(yù)報(bào)誤差可以得出正規(guī)方程 (5-4-8)或 (5-4-9)其中如果同時考慮正反向預(yù)報(bào)誤差，并對兩種預(yù)報(bào)采用同樣的自回歸系數(shù)，則可寫出含有這n個參數(shù)的2(N-n)個方程式 (5-4-

21、10)其中而實(shí)現(xiàn)的正規(guī)方程為 (5-4-11)其中 實(shí)際確定參數(shù)的過程包含兩大步驟，一是正規(guī)方程的列寫，二是正規(guī)方程的求解。一般說，對于含n個未知參數(shù)的M個超定方程組，其正規(guī)議程的建立需要次運(yùn)算，用巧列斯基（Cholesky）方法求解需要次運(yùn)算（一次運(yùn)算是指一個乘法或除法加上一個加法或減法，在比較運(yùn)算量時只考慮M和n的最高次冪）。對應(yīng)于式（5-4-10），M=2(N-n)，但在求解AR模型參數(shù)一特定情況下，由n階建立n+1階正規(guī)方程可以采用遞推的方法，下面將表明：對于單向（正向或反向）預(yù)報(bào)來說，建立正規(guī)議程的計(jì)算量為，而對于雙向預(yù)報(bào)的式（5-4-10），其計(jì)算量為。由于階數(shù)往往要從1開始搜索，

22、在合適階數(shù)未知的情況下，從1到某個階數(shù)n2逐個建立正規(guī)方程的總計(jì)算量為：以正向預(yù)報(bào)中由n=2推出n=3的情況為例，考慮到正規(guī)方程矩陣的對稱性，只寫它的下三角部分注意到中由構(gòu)成的2×2子矩陣的各元素等于中的元素減去元素所定義的內(nèi)積中的第一項(xiàng)，即最后一行中除了第一個元素外，都可由的最后一行推得，即此外，除了最后一個元素外，亦可由推得可見，除了要算一個內(nèi)積中的其余元素只要一次運(yùn)算便可求得。一般情況下，由需要多于一次的運(yùn)算。關(guān)于三種方式下正規(guī)方程的建立過程用編程格式列出如下：（1）正向預(yù)報(bào)方式算法（2）反向預(yù)報(bào)方式算法（3）雙向預(yù)報(bào)方式算法關(guān)于n×n正規(guī)方程的求解，第一步是將R（它是對稱正定陣）分解成（L為下三角陣），所需運(yùn)算量為，第二步是解，然后由，需要的運(yùn)算量為，兩步共需的運(yùn)算，這樣從到某一階數(shù)n獨(dú)立求解所需的總計(jì)算量為。由于計(jì)算中采用了三角陣的分解（巧列斯基分解），故簡稱LUD算法（Lower and Upper Triangular Matrix Decomposition）。總的說來，在情況下，全部計(jì)算中占主要的計(jì)算是建立正規(guī)方程，雖然LS-LUD算法的計(jì)算量比伯格算法大，但相差不是數(shù)量