建模與估計(jì)下

1、建模與估計(jì)講義（二） 主講人：王欣,第三節(jié) 新息序列 (Innovation sequence,定義:設(shè) (相關(guān)隨機(jī)序列正交化) 有二階矩(數(shù)學(xué)期望、方差都）的隨機(jī)序列，定義它的新息序列為： 定義: 基于k-1以前信息對(duì)y(k)的線性最小方差估計(jì) 則 規(guī)定: 定理1：新息序列是正交序列（白噪聲） 證： 射影具有無(wú)偏性 ij 不妨設(shè)ij 根據(jù)射影性質(zhì),定理2,證明: 顯然,注意,事實(shí)上,重要意義:新息序列 與原序列y(k)含有相同的統(tǒng)計(jì)信息,新息實(shí)質(zhì)上是攝影誤差,定理3: 遞推射影公式 證明,考慮線性離散定常隨機(jī)系統(tǒng) 其中狀態(tài)x(t) ,觀測(cè)y(t) ,w(t) 輸入噪聲.觀測(cè)噪聲v(t) 假設(shè)

2、1: w(t)和v(t)是零均值,方差陣各為Q和R的獨(dú)立白噪聲 即,第四節(jié) kalman濾波,假設(shè)2: x(0)與w(t)和v(t)相互獨(dú)立 Kalman濾波問(wèn)題是基于觀測(cè) (y(1)y(t) 求狀態(tài) x(j)的線性最小方差估計(jì) 極小化 即求 濾波 平滑 預(yù)報(bào),分析:有遞推射影公式 k(t+1)稱為濾波增益 對(duì)取射影 因?yàn)?事實(shí)上,新息: 引入定義,即 得 即,即 注意,簡(jiǎn)化: 即: 總結(jié):考慮隨機(jī)系統(tǒng) 問(wèn)題:基于觀測(cè) 求,Kalman Filter 初值 協(xié)方差陣,回憶:考慮隨機(jī)系統(tǒng) 問(wèn)題:基于觀測(cè) 求 定理:Kalman Filter,第四節(jié),初值: 協(xié)方差陣 例1:已知 及 求 解: 注

3、1:Kalman濾波的和預(yù)報(bào)器的不同形式 封閉形式,濾波器 預(yù)報(bào)增益 預(yù)報(bào)器,注2: Riccati方程 即: 初值: 注3:時(shí)變系統(tǒng)Kalman方程組同樣適用,有Kalman方程組 初值: 注4:Kalman濾波與RLS關(guān)系(RLS是K的特例) 考慮AR(N,可寫(xiě)成狀態(tài)空間模型 應(yīng)用Kalman Filter 定義： 具有RLS公式,注5：穩(wěn)態(tài)Kalman濾波和預(yù)報(bào)器（steady-state kalman filter and pred！ctor） Kalman方程組缺點(diǎn)：需要在線（on-line）實(shí)時(shí)（real time）計(jì)算p 問(wèn)題：P（t+1|t）?（t） 定常系統(tǒng) 即 h p為常陣

4、 定理：設(shè)系統(tǒng)是完全可觀、完全可控或 是穩(wěn)定陣 即： H rank H =n rankp p *n-1p=n H *2 H *n-1 則任意 p(t+1|t)=它滿足Riccati方程 = - H*T(H H*H+R)*(-1)H *T+PQP*T,且有關(guān)系 k（t）=k p（t|t）=p k=H H H*T+R*-1 = p *T=PQP*T p=-KH 穩(wěn)態(tài)Kalman濾波和預(yù)報(bào)器為 x(t+1|t+1)= x(t|t)+ky(t+1) =-KH x(t+1|t)= x(t|t-1)+ y(t) = -KH = k 注釋：可證 為穩(wěn)定陣 因而 x(t|t), x(t+|t)是漸近穩(wěn)定的。

5、即可任意選取初值 x(0|0)，或（0|0） 例：考慮一維系統(tǒng) x(t+1)= x(t)+bw(t) y(t) =hx(t)+v(t) Q= *2=q R= *2=r 求穩(wěn)態(tài)kalman濾波器和預(yù)報(bào)器 x(t+1|t,x(t+1|t+1)= x(t|t)+ ky(t+1) =(1-kh) x(t+1|t) = x(t|t)-1+ y(t) =-kph = k Riccati方程 = -h（hh+r）*-1h +bqb = *2- h*2(h*2+r)*-1+b*2q = *2(h*2 +r)-*2h*2/h*2 +r+b*2q (h*2 +r)= *2 r+b*2q (h*2 +r) h*2

6、*2+ r= *2 r+b*2q h*2 +b*2qr h*2 *2+(r- *2r-b*2h*2q) -b*2qr=0 作一元二次方程，取0為GT為,第五講 ARMA時(shí)間序列預(yù)報(bào)（ARMA Time series prediction) Forcasting 1、引言 預(yù)報(bào)問(wèn)題：氣象、水文、經(jīng)濟(jì)系統(tǒng)、控制,最優(yōu)預(yù)報(bào)，穩(wěn)態(tài)線性最小方差預(yù)報(bào) 穩(wěn)態(tài)： = ，已知無(wú)窮的觀測(cè)歷史 線性預(yù)報(bào)：y(t+k|t) L(y(t)y(t-1) ) 是以前歷史的線性組合 最小方差：minJ=E(y(t+k)-y(t+k|t)*2 y(t+k|t,Box-Jenkins遞推預(yù)報(bào)器 . 考慮平穩(wěn)可逆的ARMA過(guò)程y(

7、t) . A(q*1 )y(t)=c(q*1 )e(t) 其中e(t)是白噪聲：E e(t)=0 Ee(t)e(s)= *2 A(q*1 )=1+ q*1 + + q* C(q*1 )=1+ q*1 + + q* 已知（y(t)y(t-1) ）求y(t+k|t) 分析： . 平穩(wěn)性：y(t)=c(q*1 ) / A(q*1 )e(t)= e(t-j) 可逆性：e(t)=A(q*1 ) / C(q*1 )y(t)= y(t-j) 故： L(y(t)y(t-1) )=L (e(t)e(t-1) ) y(t+k|t)=proj(y(t+k)|y(t)y(t-1) ) =proj(y(t+k)|e(t

8、)e(t-1),歷史：1 wiener-kdmogorov 預(yù)報(bào)方法（1940、ARMA(p、q)=MACLO） 2 Box-Jenkins 遞推映射新方法（1970年）射影 . 3 Astrom預(yù)報(bào)方法 . 4 kalman預(yù)報(bào)方法（1960年）,一步預(yù)報(bào)： . y(t+1)=- y(t)- y(t-1) - y(t+1- ) . +e(t+1)+ e(t)+ + e(t+1- ) y(t+1|t)=- y(t)- y(t-1) - y(t+1- ) . + e(t)+ - e(t+1- ) 兩式相減得 y(t+1)- y(t+1|t)= e(t+1)(即對(duì)于平穩(wěn)ARMA過(guò)程,e(t)為y(

9、t)的新息,定理: Box-Jenkins逆推預(yù)報(bào)器 . K級(jí)預(yù)報(bào) . y(t+k)=- y(i+k-i) + e(t+k-i) ( =1) 1 k y(t+k|t)=- y(t+k-i|t) + e(t+k-i) k y(t+k|t)=- y(t+k-i|t) . 規(guī)定: y(t+k-i|t)=y(t+k-i) (t+k-i t,例1 (1-aq *1 ) y(t)=e(t) (AR(1) |a| 1 . 求 y(t+k|t)=? . 解：y(t)=ay(t-1)+e(t) . y(t+1)=ay(t)+e(t+1) . y(t|n)=ay(t) . y(t+k)=ay(t+k-1)+e(t

10、+k) (k 1) . y(t+k|t)=ay(t+k-1|t)=a*k y(t) . 例2 y(t)=(1+(q* 1 )e(t) |c| 1 . 解：y(t+1)=e(t+1)+ce(t) . y(t+1|t)=ce(t)= y(t) . y(t+1|t)=0 (k1) . y(t+1|t)+c y(t|t-1)=cy(t) . 初值：y( | 1),例3：ARMA（1.1) 解：(1aq*-1)y(t)=( 1cq-1)e(t) 求y(t+k|t)=? y(t+1)=ay(t)+e(t+1)+ce(t) y(t+1|t)=ay(t)+ ce(t) y(t+k|t)=ay(t+k-1|t)

11、 (k1) 非遞推：y(t+k|t)=a*k-1 (ay(t)+ce(t) e(t)= y(t) y(t+k|t)= a*k-1ay(t)+ y(t) = a*(k-1) y(t) = y(t,注6：穩(wěn)態(tài)Kalman濾波器寫(xiě)成Wienet的濾波器 x(t+1|t+1)=-kH x(t|t)+ky(t+1) =-kH x(t+1|t+1)= x(t|t)+ ky(t+1) x(t+1|t+1)= q-1 x(t+1|t+1) + ky(t+1) -q-1 x(t+1|t+1)= ky(t+1) 傳遞少數(shù)形式：x(t|t)= -q-1 ky(t) 例1: x(t+1)=0.5x(t)+w(t) (

12、1) Y(t)=x(t)+v(t) (2) 求：y(t+2|t) 解1：y(t+2)=x(t+2)+ v(t+2) y(t+2|t)= x(t+2|t)= 2 x(t|t) 2-0.25-1=0 1=1.1328(1.1328+1) -1 2=-0.8828舍 k=H(HHT+R)1=1.1328(1.1328+1) -1 =0.5311 K 0.5311 y(t+2|t)=0.25-y(t)=0.25-y(t) 1- q-1 1- q-110.5311 0.5 0.13275 =-y(t) 10.2344 q-1,解2：(10.5 q*-1)y(t)= w(t1)+ (10.5 q*-1)

13、v(t) (10.5 q*-1)y(t)= w(t1)+ v(t)-0.5v(t1,1+1+0.25=(1+d12) *2 d*2+4.5d+1=0 d=-0.2344 -0.5= d1 2 或d=-4.2656(舍) a(a+c) 0.5(0.5-0.2344) y(t+2/t)-y(t)= -y(t) 1+c q*-1 10.2344 q*-1 0.1328 =- y(t) 10.2344 q*-1,trm預(yù)報(bào)器 例1：一個(gè)啟發(fā)性的例子 ARMA（1.1） （1aq-1）y(t)=( 1cq-1)e(t) a 1 c 1 求y(t+2/t)=? 分析： 1+ cq*-1 y(t+2)=-e

14、(t+2) 1-aq*-1 除法綜合： 1+(c+a) q*-1,1-aq*-1 1+cq*-1 1-aq*-1,c+a) q*- 1 (c+a)q*-1a(c+a) q*-2,a(c+a) q*-2 a(c+a,y(t+2)=( 1+(c+a) q*-1-) e(t+2) 1-aq*-1,a(c+a) = e(t+2)+ (c+a)e(t+1)+ - e(t,1-aq*-1,未知 已知,a(c+a) y(t+2|t)= - e(t) 1-aq*-1 a(c+a) 1-aq*-1 = - . - y(t) 1-aq*-1 1+cq*-1 y(t+2|t)= y(t+2)- y(t+2|t)=e

15、(t+2)+(a+c) e(t+1) Ey*2(t+2|t)=1+（a+c）*2 *2,一般情況：A(q*-1) y(t)=c(q*-1) e(t) (平穩(wěn)、可逆) 求：y(t+k|t)=? c(q*-1) 分析：y(t+k)= - e(t+k) A(q*-1) c(q*-1) G(q*-1) 綜合除法: - =F(q-1) +q*-k- A(q*-1) . F(q*-1)= + q*-1+ q*-（k-1） G(q*-1)= + q*-1+ q*-ng,Diophantine方程： c(q*-1)= A(q*-1) F(q*-1)+ q*-k G(q*-1) =max（na-1,nc-k）

16、兩邊比較系數(shù)可求得F，G G(q*-1) y(t+k)= F(q*-1) e(t+k)+ - e(t) A(q*-1) G(q*-1) trm預(yù)報(bào)器 y(t+k|t)= - y(t) C(q*-1) y(t+k|t)= y(t+k)- y(t+k|t)= F(q*-1) e(t+k) Ey*2(t+k|t)= *2 *2 遞推trm預(yù)報(bào)器：C(q*-1) y(t+k|t)= G(q*-1) y(t,例2: ARMA(1,1) (1-aq*-1)y(t)=(1+cq*-1)e(t) |ak| |ck| 求 y(t+k|t)=? 解: F(q*-1) = + q*-1+ q*(k-1) =max(

17、1-1,1-k)=0 G(q*-1)= (1+cq*-1)=(1+aq*-1)( + q*-1+ q*-(k-1)+(a*-k) 比較系數(shù)法：1= =1 c= -a =c+q 0= - a =a(c+a) 0= -a =a*(k-2)(a+c) 0= -a =a*(k-1)(a+c) y(t+k|t)= y(t) 例平穩(wěn)可逆ARMA(1,2) (1-aq*-1)y(t)=(1+ q*-1+ q*-2)e(t) 求一步trmy(t+1|t)=? 及一步Box-Jenkins y(t+1|t)=,解：k=1 F(Q*-1)= =max( -k, -1)=max(2-1,1-1)=1 G(q*-1)

18、= + q*-1 Diophantine方程 1+ q*-1+ q*-2=(1-aq*-1) +q*-1( + q*-1) 1= =1 = -a + = +a = = y(t+1|t)= y(t) Box-Jenkins方法: y(t+1|t)=ay(t)+ e(t)+ e(t-1) e(t)= y(t) e(t-1)= y(t-1) = y(t) y(t+1|t)= y(t,y(t) 二者等價(jià) 封閉形式穩(wěn)態(tài)Kalman濾波器 x(t+1|t+1)=In-KH x(t|t)+ky(t) 傳遞函數(shù)形式 x(t|t)=In- q*-1*(-1)ky(t) 封閉形式穩(wěn)態(tài)Kalman預(yù)報(bào)器 x(t+1

19、|t)= In-KH x(t|t-1)+ y(t) 傳遞函數(shù)形式 x(t+1|t)=In- q*-1 y(t) 新息形式下穩(wěn)態(tài)Kalman預(yù)極器 x(t+1|t)= x(t|t) x(t|t)=x(t|t-1)+k x(t+1|t)= x(t|t-1)+k = x(t|t-1)+ 傳遞函數(shù)形式 x(t+1|t)=In-q*(-1) *(-1,第十五講 總復(fù)習(xí),基本概念 時(shí)間序列 時(shí)間序列預(yù)極 樣本函數(shù)(實(shí)現(xiàn)) Box-Jenkins 時(shí)間序列分析方法:一個(gè)樣本函數(shù)(一個(gè)實(shí)現(xiàn)) 建模 估計(jì) trm 時(shí)間序列濾波 Kalman 平穩(wěn)隨機(jī)過(guò)程: =0 =0 r(i)=E ARMA模型 (q*-1)

20、=Q(q*-1) 平穩(wěn)性 (x)=0的根在單位圓外 (q*-1)=1- q*-1+- q*-p 可逆性 (x)=0的根在單位圓外 (q*-1)=1- q*-1- q*-p 相關(guān)函數(shù) =Ez(t)z(t-k)= | | E =0 E = 標(biāo)準(zhǔn)相關(guān)函數(shù): 方差:,ARMA過(guò)程展式 成形濾波比 (規(guī)定: =1, =0 (jq) 白壓濾波比 (規(guī)定: =1, =0 (jp) =0 (jp) 注意: =0 j =0 j ARMA(p,q)=MA( )=MA( ) 充分大 AR( )=AR( ) 充分大 也可用幾何級(jí)數(shù): (1- q*-1)z(t),相關(guān)函數(shù)的計(jì)算 .最小二乘法: y(t)- y(t-1)

21、- y(t-2)- y(t-n)= y(t)= + LS結(jié)構(gòu) RLS公式: =0 p(0)=I =10*5,RELS A(q*-1)y(t)=D(q*-1) A(q*-1)=1- q*-1- q*- D(q*-1)=1- q*-1- q*- 定義: = =y(t-1)y(t- ),- - LS結(jié)構(gòu): y(t)= + 未知,用估值來(lái)代替 =y(t-1)y(t- ) - 偶合,白噪聲估值 =y(t)- 與參數(shù)估值互偶 RLS-RELS = y(t)= y(t-j)= y(t-j) (q*-1)y(t)= 多維多重RLS .射影理論 1.射影公式 proj(x|y)=x=Ex+ (y-Ey) 2.新息序列: =y(k)-proj(y(