建模与估计下

1、建模与估计讲义（二） 主讲人：王欣,第三节 新息序列 (Innovation sequence,定义:设 (相关随机序列正交化) 有二阶矩(数学期望、方差都）的随机序列，定义它的新息序列为： 定义: 基于k-1以前信息对y(k)的线性最小方差估计 则 规定: 定理1：新息序列是正交序列（白噪声） 证： 射影具有无偏性 ij 不妨设ij 根据射影性质,定理2,证明: 显然,注意,事实上,重要意义:新息序列 与原序列y(k)含有相同的统计信息,新息实质上是摄影误差,定理3: 递推射影公式 证明,考虑线性离散定常随机系统 其中状态x(t) ,观测y(t) ,w(t) 输入噪声.观测噪声v(t) 假设

2、1: w(t)和v(t)是零均值,方差阵各为Q和R的独立白噪声 即,第四节 kalman滤波,假设2: x(0)与w(t)和v(t)相互独立 Kalman滤波问题是基于观测 (y(1)y(t) 求状态 x(j)的线性最小方差估计 极小化 即求 滤波 平滑 预报,分析:有递推射影公式 k(t+1)称为滤波增益 对取射影 因为 事实上,新息: 引入定义,即 得 即,即 注意,简化: 即: 总结:考虑随机系统 问题:基于观测 求,Kalman Filter 初值 协方差阵,回忆:考虑随机系统 问题:基于观测 求 定理:Kalman Filter,第四节,初值: 协方差阵 例1:已知 及 求 解: 注

3、1:Kalman滤波的和预报器的不同形式 封闭形式,滤波器 预报增益 预报器,注2: Riccati方程 即: 初值: 注3:时变系统Kalman方程组同样适用,有Kalman方程组 初值: 注4:Kalman滤波与RLS关系(RLS是K的特例) 考虑AR(N,可写成状态空间模型 应用Kalman Filter 定义： 具有RLS公式,注5：稳态Kalman滤波和预报器（steady-state kalman filter and pred！ctor） Kalman方程组缺点：需要在线（on-line）实时（real time）计算p 问题：P（t+1|t）?（t） 定常系统 即 h p为常阵

4、 定理：设系统是完全可观、完全可控或 是稳定阵 即： 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,且有关系 k（t）=k p（t|t）=p k=H H H*T+R*-1 = p *T=PQP*T p=-KH 稳态Kalman滤波和预报器为 x(t+1|t+1)= x(t|t)+ky(t+1) =-KH x(t+1|t)= x(t|t-1)+ y(t) = -KH = k 注释：可证 为稳定阵 因而 x(t|t), x(t+|t)是渐近稳定的。

5、即可任意选取初值 x(0|0)，或（0|0） 例：考虑一维系统 x(t+1)= x(t)+bw(t) y(t) =hx(t)+v(t) Q= *2=q R= *2=r 求稳态kalman滤波器和预报器 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时间序列预报（ARMA Time series prediction) Forcasting 1、引言 预报问题：气象、水文、经济系统、控制,最优预报，稳态线性最小方差预报 稳态： = ，已知无穷的观测历史 线性预报：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递推预报器 . 考虑平稳可逆的ARMA过程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) 分析： . 平稳性：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 预报方法（1940、ARMA(p、q)=MACLO） 2 Box-Jenkins 递推映射新方法（1970年）射影 . 3 Astrom预报方法 . 4 kalman预报方法（1960年）,一步预报： . 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)(即对于平稳ARMA过程,e(t)为y(

9、t)的新息,定理: Box-Jenkins逆推预报器 . K级预报 . 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) . 规定: 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：稳态Kalman滤波器写成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) 传递少数形式：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预报器 例1：一个启发性的例子 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) (平稳、可逆) 求：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、两边比较系数可求得F，G G(q*-1) y(t+k)= F(q*-1) e(t+k)+ - e(t) A(q*-1) G(q*-1) trm预报器 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预报器：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) 比较系数法：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) 例平稳可逆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) 二者等价 封闭形式稳态Kalman滤波器 x(t+1|t+1)=In-KH x(t|t)+ky(t) 传递函数形式 x(t|t)=In- q*-1*(-1)ky(t) 封闭形式稳态Kalman预报器 x(t+1

19、|t)= In-KH x(t|t-1)+ y(t) 传递函数形式 x(t+1|t)=In- q*-1 y(t) 新息形式下稳态Kalman预极器 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)+ 传递函数形式 x(t+1|t)=In-q*(-1) *(-1,第十五讲 总复习,基本概念 时间序列 时间序列预极 样本函数(实现) Box-Jenkins 时间序列分析方法:一个样本函数(一个实现) 建模 估计 trm 时间序列滤波 Kalman 平稳随机过程: =0 =0 r(i)=E ARMA模型 (q*-1)

20、=Q(q*-1) 平稳性 (x)=0的根在单位圆外 (q*-1)=1- q*-1+- q*-p 可逆性 (x)=0的根在单位圆外 (q*-1)=1- q*-1- q*-p 相关函数 =Ez(t)z(t-k)= | | E =0 E = 标准相关函数: 方差:,ARMA过程展式 成形滤波比 (规定: =1, =0 (jq) 白压滤波比 (规定: =1, =0 (jp) =0 (jp) 注意: =0 j =0 j ARMA(p,q)=MA( )=MA( ) 充分大 AR( )=AR( ) 充分大 也可用几何级数: (1- q*-1)z(t),相关函数的计算 .最小二乘法: y(t)- y(t-1)

21、- y(t-2)- y(t-n)= y(t)= + LS结构 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结构: y(t)= + 未知,用估值来代替 =y(t-1)y(t- ) - 偶合,白噪声估值 =y(t)- 与参数估值互偶 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(