极大似然辨识及其MTLAB实现

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

2、，定义一个似然函数(1.1)L(Vl，V2, ,VM =f(ViU)f(V2U)f(VC上式的右边是n个概率密度函数的连乘， 似然函数L是e的函数。如果L达到极大值，Vk的出现概率为最大。 因此，极大似然法的实质就是求出使l达到极大值的臼的估值e。为了A便于求8 ,对式(1.1 )等号两边取对数，则把连乘变成连加，即nln L =£ ln f (Vj 0)(1.2 )i A由于对数函数是单调递增函数，当 L取极大值时，lnL也同时取极大值。求式(1.2) 对10的偏导数，令偏导数为0,可得j ln Lce=0(1.3)解上式可得e的极大似然估计Bml。2系统参数的极大似然估计设系统的

3、差分方程为a(z") y(k) = b(z)u (k) + t(k)(2.1 )式中11_na(z ) =1 - a1z- . - anz_1_1_b(z ) = b0 b1z- . - bnz因为|£(k)是相关随机向量，故(2.1 )可写成a( z") y(k) =b (z)u (k) +c( z)a(k)(2.2 )式中c( z")敏k) = E( k)(2.3 )1 一1nc(z -)=1 +c1z-+cnz(2.4 )z( k)是均值为 0的高斯分布 白噪声序 列。多项式 a(z) , b( z)和c(z)中的系数a1，a,b0L bn，6广&

4、#39;cn和序列(k)的均方差o都是未知参数。设待估参数iT6 =a，an bobn。;、(2.5)并设y(k)的预测值为AAAAAy(k) - - a1 y (k -1)-an y (k - n) b0u(k)，-，bnu(k-n)-A ,.A ,.c1 e(k 一1) +cn e(k 一 n)( 2.6 )式中e(k _i)为预测误差；ai , bi , c；为ai , bi , ci的估值。预测误差可表示为e(k) = y(k) y(k) = y(k)一.：, ai y(k i) bjU(k i)-IL i 注i =0n %I% 1n1nx cie(k -i) =(1a1z -，anz

5、 )y(k)_(b0 b1z ，bnz )u(k)_i旦-1_2_n、(ci z +c2 z+cn z )e(k)(2.7)或者1一 .1n(I ciz ，cnz )e(k)=(1az ，anz ) y (k)-z")u(k)(2.8 )(2.9 ).1(b 0 b1 z - ' - b n因此预测误差e(k)满足关系式111c(z )e(k a (z ) y( kb(z )u (k)式中nn=a( z ) = 1 a1 z ，an zb(z ) = b0 - b1 z - . - bn z、八，A c(z ) =1 c1 z -.-cn z为 *(k)与 c(z 工),假定

6、预测误差e(k)服从均值为0的高斯分布，并设序列k(k)具有相同的方差B2。因 -112a(z )和b(z )有关，所以。是被估参数8的函数。为了书写方便,把式(2.9 )写成(2.10 )(2.11 )c(z J)e(k a(z )y (kb(z J)u(k)e(k) = y(k) ' a1 y (k -1) ,' an y(k - n) - b0u (k - 1) - b1 u(k -1) 一 -bnu (k n) c1e(k 1)-cn( k n), k = n 1, n 2,或写成e(k) =y(k) -W： a/k i)-寸 0u(k i)Ge(k i)i ±

7、;i _0i _1(2.12 )-矩阵形式式中Yn =y(n 1)y(n +2)”(n + N)_|-y(n)-y(n+1) 中N =-eNeN=YN - N 71e(n 1)e(n +2)e(n +N)(2.13 )'ailan-y(1) u(n 1) . u(1)-y(2) u(n 2)u(2)e(n)e(n 1)e(1)e(2)y(n+N_1) -y(N) u(n + N) u(N )e(n N 1) e(N)因为已假定,(k)是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为112exp-一(y m)2 ；2。(2匚)(2.14 )式中y为观测值，。2和m为y的方差和均值，

8、那么12 exp-一 e (k)2 22。(2 二一-)(2.15 )对于e(k)符合高斯噪声序列的极大似然函数为L(YN 日,。)=Le(n +1),e(n +2),，e(n +N )同=f e(n +1)切 f e(n + 2)引f e(n + N )61N-2；(2二)2exp122_2.e (n 1) e (n 2)，e (n N)21N2 (2)一一， 1 _T 一、exp(-一 eNeN)2。令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这 N个方程式写成向量(2.16 )l(Yn e,s)1. (Yn '：r)T(YN -'：七)=exp-22 &qu

9、ot;2" (2= )222(2.17 )对上式(2.17)等号两边取对数得1 TNN 21 T2 eN eN) = - ln 2二-ln ；： - 了 eN eN2-222c(2.18 )n .N2 I2ln u 一 2 £ e (k)(2.19 )2- k 丑.10,可得cln L(Yn O,ct)2CCTN22c1 n.N、e2(k) =02 - k八(2.20 )LeN j 1(k)n N一 一 2、e (k)N 2 k -n 1=JN(2.21 )In L(Yn "了)=lnIn ex|(2= 2)或写为rNNIn L (Yn B,。)= In 2n -

10、 22求ln L(Yt)对。2的偏导数，令其等于式中(2.22 )。2越小越好，因为当方差 。2最小时，e2(k)最小，即残差最小。因此希望 。2的估值取最 小2(2.23 )2=min JN因为式(2.10 )可理解为预测模型，而 e(k)可看做预测误差。因此使式(2.22)最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因 此可按J最小来求a1,a,b°,bn q,Cn的估计值。由于e(k)式参数a1,a,b0,bn,c1，cn的线性函数，因此J是这些参数的二次型函数。求使ln L(YN|e,。)最大的8 ,等价于在式(2.10 )的约束条件下求

11、8使J为最小。由于J对是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下：(1) 确定初始的 说值。对于00中的a1,., a,b0,bn可按模型e(k) =a(z )y(k)-b(z )u(k)(2.24 )用最小二乘法来求，而对于S。中的6'Cn可先假定一些值。(2) 计算预测误差给出J =1 w e2(k)2 k虹并计算了2 = 1 w e2(k)(2.26 )N 口.1(3)计算J的梯度c9和海赛矩阵2斜有fe(k)(2.27 )n -N：J一 ='、e(k

12、):u j 1T任(k) 笑(k)&(k) ce(k)庄(k)&(k)ce(k) 1=|_ cai&的0犯cci如;e(k)y(k) aiy(k -1)一 any(k - n)-b0u(k) - b1u (k -1)- -bnu(k-n)-c1e(k -1)-cne(k - n)= y(kj)_C1_C2i.'aizai；:e(k - n)-cn.:ai(2.28 )cj：'e(k - j)&(2.29 )同理可得cje(k - j)A(2.30 )(2.31 )业一e(kiq cj 创acijEM将式(2.29 )移项化简，有(2.32 ),&

13、quot;_)=* ； cj4jicj5:aj 土- aij =0tai因为e(k - j)=e(k)z J(2.33 )由e（k j）求偏导，故fe(k 一 j)；：e（k）z(2.34 )将（2.34 ）代入（2.32 ),所以ny(k -i)='、j =0cje(k 一 j)nj =0fe(k)zcj :；'ai；e(k) ncjzj =0(2.35 )所以得c(z)同理可得（2.30）和(2.31 )Ac(z )-1c(z )根据（2.36 ）构造公式c( z ) = 1 - c1 zje(k)-：aije(k)-:bije(k)_n n=y(k -i)-u (k -

14、i)-e(k -i)1 ce k - (i - j)c(z )yk -(i - j) - j = y(k -i)将其代入（2.36 ）,可得.1 c(z )fek -(i - j)% fe(k)c(z ):aiaj(2.40(2.36 )(2.37 )(2.38 )(2.39 )消除c（ z汽可得;:e(k)'e(k - ij)fe(k - i 1)：a3i(2.41 )同理可得（2.37 ）和（2.38 ）式,:e(k)2(k _i j)；e(k -i)(2.42 )一'bi;：bj;：b0.:e(k)_：e(k i j)fe(k i 1)(2.43 )-：cj式(2.29)

15、、式(2.30 )和式(2.31 )均为差分方程，这些差分方程的初始条件为0,可通过求解这些差分方程，分别求出e(k)关于a1 ：., a, b0,bn, G,cn的全部偏导数，而这些偏导数分别为y(k) , u(k)和e(k)的线性函数。下面求关于 8的二阶偏导数，即：J：e(k) | ；.e(k)三=k 三 1:- 口 _ :七Tn -N-2e( k)'、e(k)rk -n 1:(2.44 )A当臼接近于真值B时，e(k)接近于0。在这种情况下，式(2.44 )等号右边第2项接近# J于0,；可近似表示为(2.45 )，一一 ：J 则利用式(2.45 )计算£J比较简单。

16、按牛顿-拉卜森计算0的新估值 奋,有 -| i； J A ；:J01 =00 1()_ (顶2浏瞄A.重复(2)至(4)的计算步骤，经过r次迭代计算之后可得 Qr ,近一步迭代计算可得(2.47 )如果则可停止计算，否则继续迭代计算。式(2.48 )表明，当残差方差的计算误差小于0.01 %时就停止计算。这一方法即使在噪声A比较大的情况也能得到较好的估计值e。3动态模型参数极大似然辨识及其MATLAB现设动态系统的模型表示为 r iiA(z)z(k) =B(z)u(k) +e(k) 1 e(k) =D(zf(k)式中，v(k)是均值为0,方差为cr2,服从正态分布的不相关随机噪声；u(k)和z

17、(k)表示系统的输入输出变量。现给出一系统模型为z(k)-1.2z(k-1)+0.6z(k-2)=u(k-1)+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);

18、endu;v=randn(1200,1); %产生正态分布随机数V=0; %计算噪声方差for i=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)=。；b1(k)=o (3);b2(k)=o(4);d1(k)=o(5);d2(k)=o (6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs