极大似然辨识及其MTLAB实现

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

2、义一个似然函数(1.1 )L(Vl，V2,V")f(V2)上式的右边是n个概率密度函数的连乘， 似然函数L是二的函数。如果L达到极大值，Vk的出现概率为最大。 因此，极大似然法的实质就是求出使L达到极大值的的估值二。为了A便于求二，对式(1.1 )等号两边取对数，则把连乘变成连加，即nIn L 八 In f (Vi v)(1.2 )i 土由于对数函数是单调递增函数，当L取极大值时，InL也同时取极大值。求式(1.2 )对二的偏导数，令偏导数为 0，可得bnu (k n) qe(k 1)-cn( k n), k = n T, n 2,(2.11)(1.3 )解上式可得二的极大似然估计V

3、W 。2系统参数的极大似然估计设系统的差分方程为a(z") y(k) = b(z)u (k) (k)(2.1 )式中a(z°) =1 pz，亠 亠anzb(z 冷=b0 Qz-.-bnz 耳因为©(k)是相关随机向量，故(2.1 )可写成(2.2 )a(z°)y(k)=b(zJ)u(k) c(z°)；(k)式中, 1、，1n/ 小,c(z-) =1 c1z - - - cnz -(2.4 )；(k)是均值为 0的高斯分布 白噪声序 列。多项式 a(z J) , b( zJ)和c(z J)中的系数 at，a,b0,bn，5，cn和序列 s(k)的

4、均方差o都是未知参数。设待估参数- a1 anb0_ bn c'" c ,0n1n(2.5)并设y(k)的预测值为AAAA.Ay(k) _ -a! y(k _ 1) 一 -an y (k _n) bou(k)亠-亠 bnu(k_ n) AAct e(k -1)亠'亠cn e(k - n)(2.6)式中e(k _i)为预测误差；ai， bi， Ci为ai， bi， ci的估值。预测误差可表示为 n nai y (k -ibiu( i)-e(k) =y(k) _y(k)二y(k) _ 八IL -nce(k -i)i z!=(1 a1Z 丄亠-亠 anz)y(k) _(bo

5、 b1 z亠-亠 bn z)u(k)_A . A .A(c1 z 亠 C2 Z 亠亠 Cnz)e(k)(2.7)或者A 1AA 1A(1 c1 z 亠'亠 cnz)e(k)=(1 a1 z 亠'亠 anz)y(k)-丄(b o A z 亠亠 b n因此预测误差:e(k)满足关系式z)u(k)(2.8)A 1A 1Ac(z )e(k) = a (z ) y( k) -b(zA)u(k)(2.9)i -1i 兰bnu (k n) qe(k 1)-cn( k n), k = n T, n 2,(2.11)bnu (k n) qe(k 1)-cn( k n), k = n T, n 2

6、,(2.11)式中AAAa( z )=1 a1 z 亠亠 an zAAAAb(z )二6亠0乙 亠'亠bnzc(z)=1 q z-.-cn z假定预测误差e(k)服从均值为0的高斯分布，并设序列(k)具有相同的方差c2。因为2(k)与c(z丄)，a(z)和b(zx)有关，所以匚2是被估参数二的函数。为了书写方便, 把式(2.9 )写成(2.10)c(z)e(k)二a(z°)y(k)b(z")u(k)e(k) = y(k) ' a1 y (k -1)，' an y(k - n) - b0u (k - 1) - b1 u(k - 1)-bnu (k n)

7、 qe(k 1)-cn( k n), k = n T, n 2,(2.11)或写成e(k) =y(k)工 ajy(k i)du(k i)qe(k i)(帚exp -厂 "2"(2 二)2对上式(2.17 )等号两边取对数得12)i ±i _0i _1令k=n+1,n+2,n+N,可得e(k)的N个方程式，把这 N个方程式写成向量-矩阵形式eN(2.13 )式中ail2(2.17)22(2.17)2Yn =y(n 1)y(n +2)e(n 1)e(n +2)eNanboGn屮 n + N)_|-y( n)-y( n +1)y(1)u(n 1) -u(1)e(n)e(

8、1)y(2)u(n 2)u(2)e(n 1)ae(2)y(N)u(n N)u(N)e(n N-1)e( N)-1)八_-y( n 十 N因为已假定怡(k)是均值为o的高斯噪声序列，高斯噪声序列的概率密度函数为1 12?exp 2(y m)2 i2a(2栄)(2.14)式中y为观测值，：2和m为y的方差和均值，那么1 1 2 f?exp丁e (k)2 22 cr(2 二一-)对于e(k)符合高斯噪声序列的极大似然函数为(2.15)L(YN 日,<!)=Le(n +1),e(n +2),，e(n +N )|日=f e(n +1)0 f e(n +2)日f e(n + N )日1N2 &quo

9、t;2 (2匚)2expI 2222e (n 1) e (n 2)亠(n N)2 cr1N2 "2" (2 二)1 Texp( pNeN)2 CT2(2.17)22(2.17)2(2.16)L(Yn B,g1(Yn -日)T(Yn -日)、2(2.17)22In L(Yn n) =ln1N(2= 2)三1In exp(22 CTtNNeNe”)二In 2二In2 2(2.18 )e(k)二 y(k) - y(k)(2.25 )2或写为e(k)二 y(k) - y(k)(2.25 )2e(k)二 y(k) - y(k)(2.25 )2In L(Yn 日，r)NIn 2 二2I

10、nn .Nv e2(k)k史-1(2.19 )求In L(Yn)；)对匚2的偏导数，令其等于0,可得cln L(Yn 0,cr)2CCT1 n .N' e2(k) =02-k 艾 1(2.20 )(k)n N2Z e (k)N 2 k -n 1(2.21 )e(k)二 y(k) - y(k)(2.25 )2式中e(k)二 y(k) - y(k)(2.25 )2e(k)二 y(k) - y(k)(2.25 )2(2.22 )c2越小越好，因为当方差 小c 2最小时，e2(k)最小，即残差最小。因此希望:二2的估值取最2A 2 二=min JN因为式(2.10 )可理解为预测模型，而e(k

11、)可看做预测误就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。此可按J最小来求ai v., a, b0bn, Ci,差。因此使式(2.23 )2.22 )最小因Cn的估计值。e(k)二 y(k) - y(k)(2.25 )2由于e(k)式参数a1，,a,b0,bn,c1cn的线性函数，因此J是这些参数的二次型函数。AA求使In L(Yn二二)最大的二，等价于在式(2.10 )的约束条件下求 二使J为最小。由于J对C是非线性的，因而求J的极小值问题并不好解， 只能用迭代方法求解。 求J极小值的常用 迭代算法有拉格朗日乘子法和牛顿 -拉卜森法。下面介绍牛顿-拉卜森法。

12、整个迭代计算步骤 如下：(1) 确定初始的00值。对于00中的a1,a,b0,bn可按模型e(k)二 a(z °)y (k) b(z)u( k)(2.24 )用最小二乘法来求，而对于氏中的Cn可先假定一些值。计算预测误差e(k)二 y(k) - y(k)(2.25 )给出1 n NJ = _、 e2 J2 k仝1(k)并计算2CTn -Nzk -n(k)(2.26)计算J的梯度V和海赛矩阵2了有式中n -N；Je(k):二 心1：e(k)(2.27)fe(k) _ ：e(k)；:e(k)：e(k)'v a.:an.:bo.:e(k)：e(k);：GTce(k)1CCn;e(k

13、).:aiy(k) ' ay(k -1)亠亠any(k _n) _b0u(k) _du(k -1)_bnu(k _ n) -c,e(k -1)-cne(k - n)ce(k -1)y (k - i) - C1- - C2a:：e(k -2)；=e(k 一 n)-Cna(2.28)型=y(k-i)，、aijCj：e(k - j)a(2.29)同理可得e(k - j)(2.30)k)_e(k-i)£ Cjj =i;：e(k - j)(2.31)Ci将式(2.29 )移项化简，有ny(k -i)S j ±Cje(k - j)八:问j =e：e(k- j)；a iCj(2.

14、32)(2.33 )(2.34 )因为e(k _ j) = e(k) z由e(k j)求偏导，故:e(kj) ：e(k)z同理可得（2.37 ）和（2.38 ）式同理可得（2.37 ）和（2.38 ）式将（2.34 ）代入（2.32 ）,所以n：e（k 一 j）n：e（k）z：e(k) /y(k -i)=匚 j八Cjcjzj =0；：aij =0；：ai；'a ij 三c(z)=1 qzcn z(2.35 )同理可得（2.37 ）和（2.38 ）式同理可得（2.37 ）和（2.38 ）式所以得(2.36 )（x :e（k）c（z ）y（k -i）矽i同理可得（2.30 ）和（2.31

15、）为c(ze(k)巾i-u (k - i)(2.37 )c(z严&i-e(k -i)(2.38 )同理可得（2.37 ）和（2.38 ）式同理可得（2.37 ）和（2.38 ）式根据（2.36 ）构造公式(2.39 )c（z）：ek（i - j打二 丫出（i 一） j= 丫代 _i）将其代入（2.36 ），可得c(zgk -(i - j)= c(z)：e(k)'ai(2.40同理可得（2.37 ）和（2.38 ）式同理可得（2.37 ）和（2.38 ）式消除c（ z °）可得同理可得（2.37 ）和（2.38 ）式同理可得（2.37 ）和（2.38 ）式：e(k i

16、j)&j：e(k - i 1)ca1(2.41 )同理可得（2.37 ）和（2.38 ）式:e(k)：e(k -i j)：e(k -i)(2 42 )-:bicbjcbo；e(k)：e(k _ij)：e(k i 1)(2.43 )式(2.29 )、式(2.30 )和式(2.31 )均为差分方程，这些差分方程的初始条件为0,可通过求解这些差分方程，分别求出e(k)关于ai,.,a,bo,bn,G,5的全部偏导数，而这些偏导数分别为 y (k)，u (k)和 e(k)的线性函数。下面求关于二的二阶偏导数，即：e(k) |l MlTJ (k S2e(k)(k)-k -n 1:二(2.44)A

17、当二接近于真值 二时，e(k)接近于0。在这种情况下，式(2.44 )等号右边第2项接近nN :e(k) :e(k) Tk ± 1:宀 _:-二于0，、可近似表示为(2.45)则利用式(2.45 )计算厶2比较简单。按牛顿-拉卜森计算r的新估值R，有(2.46 )r2J(7Z7)重复(2)至(4)的计算步骤，经过r次迭代计算之后可得R，近一步迭代计算可得如果C0(2.47 ):10_4(2.48 )二 r则可停止计算，否则继续迭代计算。式(2.48 )表明，当残差方差的计算误差小于0.01 %时就停止计算。这一方法即使在噪声A比较大的情况也能得到较好的估计值二。3动态模型参数极大似然

18、辨识及其MATLA实现设动态系统的模型表示为广iiA(z)z(k) =B(z)u(k) +e(k)1e(k)二D(z：v(k)式中，v(k)是均值为0，方差为二，服从正态分布的不相关随机噪声；u(k)和z(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(

19、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=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*one

20、s(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