极大似然辨识及其MTLAB实现

极大似然辨识及其MATLAB实现摘要：极大似然参数估计方法是以观测值的出现概率为最大作为准则的，这是一种很普遍的参数估计方法，在系统辨识中有着广泛的应用。本文主要探讨了极大似然参数估计方法以及动态模型参数的极大似然辨识并且对其进行了MATLAB实现。关键词：极大似然辨识MATLAB仿真迭代计算1极大似然原理设有离散随机过程}{kV与未知参数有关，假定已知概率分布密度)(kVf。如果我们得到n个独立的观测值,21,VV…nV,，则可得分布密度)(1Vf，)(2Vf，…,)(nVf。要求根据这些观测值来估计未知参数，估计的准则是观测值{}{kV}的出现概率为最大。为此，定义一个似然函数)()()(),,,(2121nnVfVfVfVVVL（1.1）上式的右边是n个概率密度函数的连乘，似然函数L是的函数。如果L达到极大值，}{kV的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的的估值。为了便于求，对式（1.1）等号两边取对数，则把连乘变成连加，即niiVfL1)(lnln（1.2）由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式（1.2）对的偏导数，令偏导数为0，可得0lnL（1.3）解上式可得的极大似然估计ML。2系统参数的极大似然估计设系统的差分方程为)()()()()(11kkuzbkyza（2.1）式中111()1...nnazazaz1101()...nnbzbbzbz因为)(k是相关随机向量，故（2.1）可写成)()()()()()(111kzckuzbkyza（2.2）式中)()()(1kkzc（2.3）nnzczczc1111)(（2.4）)(k是均值为0的高斯分布白噪声序列。多项式)(1za，)(1zb和)(1zc中的系数nnccbbaa,,,,,10,1和序列)}({k的均方差都是未知参数。设待估参数naa1[nbb0Tncc1（2.5）并设)(ky的预测值为)()()()1()(01nkubkubnkyakyakynn)()1(1nkeckecn（2.6）式中)(ike为预测误差；ia，ib，ic为ia，ib，ic的估值。预测误差可表示为)()()()()()(01ikubikyakykykykeniinii)()()()1()(110111kuzbzbbkyzazaikecnnnnnii)()(2211kezczczcnn（2.7）或者)()1(11kezczcnn=)()1(11kyzazann)()(110kuzbzbbnn（2.8）因此预测误差)(ke满足关系式)()()()()()(111kuzbkyzakezc（2.9）式中nnzazaza1111)(nnzbzbbzb1101)(nnzczczc1111)(假定预测误差)(ke服从均值为0的高斯分布，并设序列)(ke具有相同的方差2。因为)(ke与)(1zc，)(1za和)(1zb有关，所以2是被估参数的函数。为了书写方便，把式（2.9）写成)()()()()()(111kuzbkyzakezc（2.10）)1()1()()1()()(101kubkubnkyakyakyken,2,1),()1()(1nnknkckecnkubnn（2.11）或写成)()()()()(101ikecikubikyakykeniiniinii（2.12）令k=n+1,n+2,…,n+N,可得)(ke的N个方程式，把这N个方程式写成向量-矩阵形式NNNYe（2.13）式中)()2()1(NnynynyYN,)()2()1(NneneneeN,nnbbaa01)1()1()(NnynynyN)()2()1(Nyyy)()2()1(Nnununu)()2()1(Nuuu)1()1()(Nnenene)()2()1(Neee因为已假定)(ke是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为])(21exp[)2(122212myf（2.14）式中y为观测值，2和m为y的方差和均值，那么)](21exp[)2(122212kef（2.15）对于)(ke符合高斯噪声序列的极大似然函数为)21exp()2(1)]}()2()1([21exp{)2(1])([])2([])1([])(,),2(),1([),(222222222NTNNNNeeNneneneNnefnefnefNneneneLYL（2.16）或]2)()(exp[)2(1),(222NTNNNYYYL（2.17）对上式（2.17）等号两边取对数得NTNNTNNNeeNNeeYL2222221ln22ln2)21exp(ln)2(1ln),(ln（2.18）或写为NnnkNkeNNYL1222)(21ln22ln2),(ln（2.19）求),(lnNYL对2的偏导数，令其等于0，可得0)(212),(ln12422NnnkNkeNYL（2.20）则JNkeNkeNNnnkNnnk2)(212)(112122（2.21）式中NnnkkeJ12)(21（2.22）2越小越好，因为当方差2最小时，)(2ke最小，即残差最小。因此希望2的估值取最小JNmin22（2.23）因为式（2.10）可理解为预测模型，而e(k)可看做预测误差。因此使式（2.22）最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因此可按J最小来求nnccbbaa,,,,,10,1的估计值。由于e(k)式参数nnccbbaa,,,,,10,1的线性函数，因此J是这些参数的二次型函数。求使),(lnNYL最大的，等价于在式（2.10）的约束条件下求使J为最小。由于J对ic是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下：(1)确定初始的0值。对于0中的nbbaa,,,0,1可按模型)()()()()(11kuzbkyzake（2.24）用最小二乘法来求，而对于0中的ncc,1可先假定一些值。(2)计算预测误差)()()(kykyke（2.25）给出NnnkkeJ12)(21并计算NnnkkeN122)(1（2.26）(3)计算J的梯度J和海赛矩阵22J,有)()(1kekeJNnnk（2.27）式中nakeakeke)()()(1nbkebke)()(0Tnckecke)()(1)()1()()()1()([)(101nkubkubkubnkyakyakyaakennii)]()1(1nkeckecniniiankecakecakeciky)()2()1()(21（2.28）即injjiajkecikyake)()()(1（2.29）同理可得injjibjkecikubke)()()(1（2.30）injjicjkecikecke)()()(1（2.31）将式（2.29）移项化简，有injjinjjiajkecajkecakeiky)()()()(01（2.32）因为jzkejke)()(（2.33）由)(jke求偏导，故ijiazkeajke)()(（2.34）将（2.34）代入（2.32），所以jnjjiijnjjinjjzcakeazkecajkeciky000)()()()(（2.35）nnzczczc1111)(所以得)()()(1ikyakezci（2.36）同理可得（2.30）和（2.31）为)()()(1ikubkezci（2.37）)()()(1ikeckezci（2.38）根据（2.36）构造公式)(])([)]([)(1ikyjjikyajikezcj（2.39）将其代入（2.36），可得ijakezcajikezc)()()]([)(11(2.40）消除)(1zc可得1)1()()(aikeajikeakeji（2.41）同理可得（2.37）和（2.38）式0)()()(bikebjikebkeji（2.42）1)1()()(cikecjikeckeji（2.43）式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为0，可通过求解这些差分方程，分别求出e(k)关于nnccbbaa,,,,,10,1的全部偏导数，而这些偏导数分别为)}({ky，)}({ku和)}({ke的线性函数。下面求关于的二阶偏导数，即NnnkNnnkkekekekeJT122122)()()()(（2.44）当接近于真值时，e(k)接近于0。在这种情况下，式（2.44）等号右边第2项接近于0，22J可近似表示为TNnnkkekeJ)()(122（2.45）则利用式（2.45）计算22J比较简单。(4)按牛顿-拉卜森计算的新估值1，有021201)(JJ（2.46）重复(2)至(4)的计算步骤，经过r次迭代计算之后可得r，近一步迭代计算可得rJJrr121)(2（2.47）如果4221210rrr（2.48）则可停止计算，否则继续迭代计算。式（2.48）表明，当残差方差的计算误差小于%01.0时就停止计算。这一方法即使在噪声比较大的情况也能得到较好的估计值。3动态模型参数极大似然辨识及其MATLAB实现设动态系统的模型表示为)()()()()()()()(111kvzDkekekuzBkzzA式中，)(kv是均值为0，方差为2，服从正态分布的不相关随机噪声；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.2v(k-2)其中v(k)为随机信号，输入信号是幅值为1的M系列或随机信号，试用递推的极大似然法求系统辨识的参数Φ。程序如下