极大似然辨识及其MTLAB实现

1、极大似然辨识及其MATLAB实现摘 要：极大似然参数估计方法是以观测值的出现概率为最大作为准则的，这是一种很普遍 的参数估计方法，在系统辨识中有着广泛的应用。本文主要探讨了极大似然参数估计方法以 及动态模型参数的极大似然辨识并且对其进行7MATLAB实现。关键词：极大似然辨识MATLAB仿真迭代计算1极大似然原理设有离散随机过程匕与未知参数9有关，假定已知概率分布密度f (p)。如果我们 得到n个独立的观测值V ,V,V，则可得分布密度f (V9 )，f (V 9)，f (V 9)。12, n12n要求根据这些观测值来估计未知参数9，估计的准则是观测值V 的出现概率为最大。k为此，定义一个似然

2、函数l(vv-，v)= f (匕 9)f (V29)-f (V)()上式的右边是n个概率密度函数的连乘，似然函数L是9的函数。如果L达到极大值，匕 的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的9的估值9。为了 便于求9，对式(1.1)等号两边取对数，则把连乘变成连加，即(1.2)ln L = ln f (V 9 )(1.2)i i =1由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式(1.2)对9的偏导数，令偏导数为0，可得公=0(1.3)89(1.3)解上式可得9的极大似然估计ML。2系统参数的极大似然估计设系统的差分方程为a(z-1)y(k) =

3、b(z-1)u(k) + & (k)(2.1)式中a(z-1) = 1 + a z-1 +. + a z-nb(z-1) = b + b z-1 +. + b z-n因为h (k)|是相关随机向量，故(2.1)可写成a(z-1)y(k) = b(z-1)u(k) + c(z-1)e (k)(2.2)式中(2.6)c的估值。预测误差可表示为e(2.6)c的估值。预测误差可表示为e(k ) = y (k)-y (k ) = y (k)-&y (k - i) + &u (k - i) +i=1i=0(2.4)c ( z-1) = 1 + c z T HF c z -n(2.4)e (k)是均值为0的

4、高斯分布白噪声序列。多项式 a(z-1) , b(z-1)和c(z-1)中的系数 a ,a,b， b ,c，c和序列e (k)的均方差n都是未知参数。1,0 n 1 n设待估参数 TOC o 1-5 h z 。=a a b b c c 】T(2.5)1 n 0 n 1 n并设y(k)的预测值为y (k) = - a y (k -1)a y (k - n) + b u (k) + b u (k - n) +1nc e(k -1) + c e(k - n)1n式中e(k - i)为预测误差；a.，b.，c.为a,&e(k-i&e(k-i)ii=1=(1 + a z-1 + an z-n)y(k)

5、- (b0+b z-1 + bn z-n)u(k)- TOC o 1-5 h z (g z-1 + c2 z-2 + + cn z-n)e(k)(2.7)或者(1 + 匕 z-1 + cz-n )e(k) = (1 + a z-1 + a z-n) y(k)-(b0 + b z-1 + + bn z-n)u(k)(2.8)因此预测误差b(k异满足关系式c (z T)e(k) = a (z-1) y (k) 一 b (z -1)u (k)(2.9)式中A / 一 -1.A一.Aa(z-1) = 1 + az-1 + az-nb (z -1) = b + bz -1 + bz-nA / 一 -1.

6、A一.Ac(z-1) = 1 + c z-1 + c z-n假定预测误差e(k)服从均值为0的高斯分布，并设序列b(k)具有相同的方差6。因 为e(k)与c(z-1)，a(z-1)和b(z-1)有关，所以g是被估参数。的函数。为了书写方便， 把式(2.9)写成c (z-1)e( k) = a (z-1) y (k) - b( z-1)u (k)(2.10)e(k) = y(k) + a y(k -1) + a y(k -n) -b u(k -1) -bu(k -1)1n01b u(k - n) - c e(k -1)c (k - n),k = n +1,n + 2,.(2.11)n1n或写成e

7、(k) = y( k) +a y(k - i)-ib u (k i)ice(k i)i(2.12)i=1i=0ie(k) = y( k) +a y(k - i)-ib u (k i)ice(k i)i(2.12)i=1i=0i =1令k=n+1,n+2,n+N,可得e（k）的N个方程式，把这N个方程式写成向量-矩阵形式（2.13）式中y (n +1)e(n +1)a -1:Y =Ny (n + 2) :,e(n + 2):,0 =ab:0_ y (n + N)_e(n + N)b _N=中 n9y (n)y (n +1)y(1)y(2)u (n +1)u (1) u(n + 2) u(2)e(

8、n)e( n +1)-e(1)e(2)u (n + N)u (N)e(n + N 一 1)e(N)y （n + N 1） - y （N）因为已假定异是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为f =1exp -（y -m）22a2（2 兀b 2）2（2.14）式中y为观测值，b 2和m为y的方差和均值，那么11f =exp e2（k）、12a 2（2兀a 2）2对于e（k）符合高斯噪声序列的极大似然函数为（2.15）L(Yn P，a) = Le(n +1), e(n + 2),e(n + N) p = f e( n +1) 0 f e(n + 2) Q f e(n + N) Q ex

9、p-e2(n +1) + e 2(n + 2) + e2(n + N) = , n2a 2n(2 兀a 2)2(2 兀a 2)2exp(eTe ) 2a 2 n n或1(Y O0)t(Y 20)】L(Y 0 ,a) =exp n n(2兀a 2)2（2.16）（2.17）)一N)一Nln2K-Nlnb2 - 1N 22er e2b 2 N N(2.18)(2.19)(2.20)(2.21)(2.22)1、，1InL(Y p,b) = In+ Inexp(一ere(2 兀b 2)2或写为NN1 n+NlnL(Y p,b) = ln2兀一Inb2 -罗 e2(k) k =n+1求ln(0,b)对b

10、2的偏导数，令其等于0,可得6ln心，b)= + 上罗Ne2(k) = 0 db 22b 2 2b 4k=n+1则2=4 勿e2(k) = 1 勿e2(k) = &NN 2Nk =n+1k =n+1式中J = N e 2(k)2k=n+1b2越小越好，因为当方差b2最小时，e2(k)最小，即残差最小。因此希望b 2的估值取最(2.23)=min J N(2.23)因为式(2.10)可理解为预测模型，而e(k)可看做预测误差。因此使式(2.22)最小 就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因 此可按J最小来求a ,a,b，b ,c，c的估计值。 1,0 n

11、 1 n由于e(k)式参数a ,a,b， b ,c，c的线性函数，因此J是这些参数的二次型函数。 1,0 n 1 n求使lnL(Y 0,b)最大的(A，等价于在式(2.10)的约束条件下求谷使J为最小。由于J对 Nc.是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤 如下：(1)确定初始的00值。对于00中的,a,b0,bn可按模型e(k) = a (z -1) y (k) - b (z-1)u (k)(2.24)用最小二乘法来求，而对于P0中的c1,cn可先假定一些值。(2)计算预

12、测误差(2.25)e( k) = y (k) - y (k)(2.25)给出1寸NJ - e2（k）2k =n+1并计算2 -艺e2（k）Nk = n+1dJFd 2 J（3）计算J的梯度所 和海赛矩阵福厂，有 囤2(2.26)式中de(k) P de(k)=60da1-1de(k)dadadJ n+ n 、合e(k)=乙 e(k) d0d0（2.27）k=n+16e(k) de(k) de(k) 6e(k) dadbdb6cde(k)dcny(k) + a y(k -1) + a y(k -n) -b u(k) -bu(k -1)b u(k -n)-c e(k -1)c e(k - n)（2

13、.28）de(k -1) de(k - 2)de(k - n)（2.28）=y(I)-七 -七-CnTiii即同理可得de(k)dai=y(k - i)-乙、jdaj=1i(2.29)de(k)了de(k - j)=-u (k - i) 一 c db. 1 j db（2.30）住= -e(k- i)-乙巡立dcjdcij=1i（2.31）将式（2.29）移项化简，有y(k - i)=岑 + lLc= 七ij=1ij=i(2.32)因为e(k j) = e(k )z - j由e(k j)求偏导，故（2.33）de(k j) _ de(k)zjdadaii将(2.34)代入(2.32)，所以（2.

14、34）y(k乙 de)二二 de=些, jdajdadajj=oij=0ii j=0（2.35）c( z-1) = 1 + C1 z1 + c z n所以得c( z-1)de(k = y (k i) dai同理可得(2.30)和(2.31)为（2.36）de(k)c( z i)= u (k i)dbi(.de( k)c( zi)= e(k i)dci根据(2.36)构造公式c(z-1)de伙(i j) = yk (i j) j = y(k i) daj将其代入(2.36)，可得c(z-1)dek -(i - j) = c (z-1)业dadaji（2.37）（2.38）（2.39）(2.40)

15、消除c（z-1）可得（2.41）dadadaij1de(k)de(k i + j)de( k i +（2.41）dadadaij13e(k)_ 3e(k - i + j)_ 3e(k - i)（2.42）3b3b3bij03e(k)3e(k 一 i + j)3e( k 一 i +1)（2.43）3ci3cj3c1可通过求解这些差分方程，式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为0,可通过求解这些差分方程，分别求出e（k）关于a ,a,b， b ,c，c的全部偏导数，而这 1,0 n 1 n些偏导数分别为y（k），u（k）和e（k）的线性函数。下面求关于

16、。的二阶偏导数，即T禁=留警催?+留仆潘（2.44）k=n+1k=n+1当接近于真值。时，e（k）接近于0。在这种情况下，式（2.44）等号右边第2项接近于0，可近似表示为于0，可近似表示为（2.45）d2J _ 艺n 8e(k)8e(k) T（2.45）30 230 |_ 30 _|k _n+1则利用式（2.45）计算则利用式（2.45）计算32 J30 2比较简单。（4）按牛顿-拉卜森计算0的新估值# 1，有01 _001 _0 0 -(邑)-1 -J30 230（2.46）00重复（2）至（4）的计算步骤，经过r次迭代计算之后可得&，近一步迭代计算可得0 =0 0 =0 -r+1r(皂)

17、-1-丑30 /30（2.47）如果品-G品-G2r+1r 10 - 4(2.48)G2 r则可停止计算，否则继续迭代计算。式(2.48)表明，当残差方差的计算误差小于.01%时就停止计算。这一方法即使在噪声 比较大的情况也能得到较好的估计值#。3动态模型参数极大似然辨识及其MATLAB实现设动态系统的模型表示为J A( z-1) z (k) = B( z-1) (k) + e(k)e(k) = D (z -1)v(k)式中，v(k)是均值为0，方差为。2服从正态分布的不相关随机噪声；u(k)和z(k)表示系统的输入输出变量。现给出一系统模型为z(k)T.2z(kT)+0.6z(k-2)=u(

18、kT)+0.5(k-2)+e(k)e(k)=v(k)- v(k-1)+0.2 v(k-2)其中v(k)为随机信号，输入信号是幅值为1的M系列或随机信号，试用递推的极大似然法求系统辨识的参数中。程序如下：cleara(1)=1;b(1)=0;d(1)=0;u(1)= d(1);z(1)=0;z(2)=0;%初 始化for i=2:1200%产生 m 序列 u(i)a(i)=xor(c(i-1),d(i-1);b(i)=a(i-i);c(i)=b(i-1);d(i)=c(i-1);u(i)=d(i);endu;v=randn(1200,1);%产生正态分布随机数V=0; %计算噪声方差for i=

19、1:1200V=V+v(i)*v(i);endV1=V/1200;for k=3:1200 %根据v和u计算zz(k)=1.2*z(k-1)-0.6*z(k-2)+u (k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);endo1=0.001*ones(6,1);p0=eye(6,6);卿幅初值zf(1)=0.1; zf(2)=0.1; vf(2)=0.1; vf(1)=0.1; uf(2)=0.1; uf(1)=0.1;%迭代计算参数值和误差值for k=3:1200h=-z(k-1); -z(k-2); u(k-1); u(k-2); v(k-1); v(k-2);hf=h;K=p0*hf*inv(hf*p0*hf+1);p=eye(6,6)-K*hf*p0;v(k)= z(k)-h*o1;o=o1+K*v(k);p0=p;o1=o;a2(k)=o(2);b1(k)=o(3);b2(k)=o(4);d1(k)=o(5) ;d2(k)=o(6) ;e1(k)