极大似然辨识及其MTLAB实现

极大似然辨识及其MATLAB 实现摘 要：极大似然参数估计方法是以观测值的出现概率为最大作为准则的，这是一种很普遍的参数估计方法，在系统辨识中有着广泛的应用。

本文主要探讨了极大似然参数估计方法以及动态模型参数的极大似然辨识并且对其进行了MATLAB 实现。

关键词：极大似然辨识 MATLAB 仿真 迭代计算1 极大似然原理设有离散随机过程}{k V 与未知参数θ有关，假定已知概率分布密度)(θk V f 。

如果我们得到n 个独立的观测值,21,V V …n V ,，则可得分布密度)(1θV f ，)(2θV f ，…,)(θn V f 。

要求根据这些观测值来估计未知参数θ，估计的准则是观测值{}{k V }的出现概率为最大。

为此，定义一个似然函数)()()(),,,(2121θθθθn n V f V f V f V V V L = （1.1）上式的右边是n 个概率密度函数的连乘，似然函数L 是θ的函数。

如果L 达到极大值，}{k V 的出现概率为最大。

因此，极大似然法的实质就是求出使L 达到极大值的θ的估值∧θ。

为了便于求∧θ，对式（1.1）等号两边取对数，则把连乘变成连加，即 ∑==ni iV f L 1)(ln ln θ （1.2）由于对数函数是单调递增函数，当L 取极大值时，lnL 也同时取极大值。

求式（1.2）对θ的偏导数，令偏导数为0，可得ln =∂∂θL（1.3）解上式可得θ的极大似然估计ML ∧θ。

2 系统参数的极大似然估计设系统的差分方程为)()()()()(11k k u z b k y z a ξ+=-- （2.1） 式中111()1...nn a z a z a z ---=+++1101()...nn b z b b z b z ---=+++因为)(k ξ是相关随机向量，故（2.1）可写成)()()()()()(111k z c k u z b k y z a ε---+= （2.2） 式中)()()(1k k z c ξε=- （2.3）nn z c z c z c ---+++= 1111)( （2.4）)(k ε是均值为0的高斯分布白噪声序列。

多项式)(1-z a ，)(1-z b 和)(1-z c 中的系数n n c c b b a a ,,,,,10,1和序列)}({k ε的均方差σ都是未知参数。

设待估参数n a a 1[=θ n b b 0 ]Tn c c 1 （2.5） 并设)(k y 的预测值为+-+++-----=∧∧∧∧∧)()()()1()(01n k u b k u b n k y a k y a k y n n)()1(1n k e c k e c n -++-∧∧（2.6） 式中)(i k e -为预测误差；i a ∧，i b ∧，i c ∧为i a ，i b ，i c 的估值。

预测误差可表示为+-+-⎢⎣⎡--=-=∑∑=∧=∧∧)()()()()()(01i k u b i k y a k y k y k y k e n i i n i i-+++-+++=⎥⎦⎤--∧-∧∧-∧-∧=∧∑)()()()1()(110111k u z b z b b k y z a z a i k e c nn n n ni i )()(2211k e z c z c z c n n -∧-∧-∧+++ （2.7）或者)()1(11k e z c z c n n -∧-∧+++ =-+++-∧-∧)()1(11k y z a z a nn)()(110k u z b z b b nn -∧-∧∧+++ （2.8） 因此预测误差{})(k e 满足关系式)()()()()()(111k u z b k y z a k e z c -∧-∧-∧-= （2.9） 式中n n z a z a z a -∧-∧-∧+++= 1111)( n n z b z b b z b -∧-∧∧-∧+++= 1101)( n n z c z c z c -∧-∧-∧+++= 1111)(假定预测误差)(k e 服从均值为0的高斯分布，并设序列{})(k e 具有相同的方差2σ。

因为{})(k e 与)(1-∧z c ，)(1-∧z a 和)(1-∧z b 有关，所以2σ是被估参数θ的函数。

为了书写方便，把式（2.9）写成)()()()()()(111k u z b k y z a k e z c ----= （2.10）-------++-+= )1()1()()1()()(101k u b k u b n k y a k y a k y k e n,2,1),()1()(1++=------n n k n k c k e c n k u b n n （2.11）或写成)()()()()(11i k e c i k u b i k y a k y k e ni in i in i i-----+=∑∑∑=== （2.12）令k=n+1,n+2,…,n+N,可得)(k e 的N 个方程式，把这N 个方程式写成向量-矩阵形式θN N N Y e Φ-= （2.13） 式中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++=)()2()1(N n y n y n y Y N ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++=)()2()1(N n e n e n e e N ,⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=n n b b a a 01θ⎢⎢⎢⎢⎣⎡-+-+--=Φ)1()1()(N n y n y n y N )()2()1(N y y y --- )()2()1(N n u n u n u +++ )()2()1(N u u u)1()1()(-++N n e n e n e ⎥⎥⎥⎥⎦⎤)()2()1(N e e e因为已假定{})(k e 是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为])(21ex p[)2(122212m y f --=σπσ （2.14）式中y 为观测值，2σ和m 为y 的方差和均值，那么)](21ex p[)2(122212k e f σπσ-=（2.15） 对于)(k e 符合高斯噪声序列的极大似然函数为)21exp()2(1)]}()2()1([21exp{)2(1])([])2([])1([])(,),2(),1([),(222222222NT N NNN e e N n e n e n e N n e f n e f n e f N n e n e n e L Y L σπσσπσθθθθσθ-=++++++-=+++=+++=（2.16）或]2)()(exp[)2(1),(222σθθπσσθΦ-Φ--=N T N NN Y Y Y L （2.17）对上式（2.17）等号两边取对数得N T N NT N N N e e N N e e Y L 2222221ln 22ln 2)21ex p(ln )2(1ln),(ln σσπσπσσθ---=-+= （2.18）或写为∑++=---=N n n k N k e N N Y L 1222)(21ln 22ln 2),(ln σσπσθ （2.19） 求),(ln σθN Y L 对2σ的偏导数，令其等于0，可得0)(212),(ln 12422=+-=∂∂∑++=Nn n k N k e N Y L σσσσθ （2.20）则J N k e N k e NN n n k Nn n k 2)(212)(112122===∑∑++=++=∧σ （2.21） 式中∑++==N n n k k e J 12)(21 （2.22）2σ越小越好，因为当方差2σ最小时，)(2k e 最小，即残差最小。

因此希望2σ的估值取最小J Nmin 22=∧σ （2.23） 因为式（2.10）可理解为预测模型，而e(k)可看做预测误差。

因此使式（2.22）最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。

因此可按J 最小来求n n c c b b a a ,,,,,10,1的估计值。

由于e(k)式参数n n c c b b a a ,,,,,10,1的线性函数，因此J 是这些参数的二次型函数。

求使),(ln σθN Y L 最大的∧θ，等价于在式（2.10）的约束条件下求∧θ使J 为最小。

由于J 对i c 是非线性的，因而求J 的极小值问题并不好解，只能用迭代方法求解。

求J 极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。

下面介绍牛顿-拉卜森法。

整个迭代计算步骤如下：(1)确定初始的0∧θ值。

对于0∧θ中的n b b a a ,,,0,1可按模型)()()()()(11k u z b k y z a k e -∧-∧-= （2.24） 用最小二乘法来求，而对于0∧θ中的nc c ,1可先假定一些值。

(2)计算预测误差)()()(k y k y k e ∧-= （2.25）给出∑++==N n n k k e J 12)(21并计算∑++=∧=Nn n k k eN 122)(1σ （2.26）(3)计算J 的梯度θ∂∂J和海赛矩阵 22θ∂∂J,有θθ∂∂=∂∂∑++=)()(1k e k e J N n n k （2.27） 式中⎢⎣⎡∂∂∂∂=∂∂n a k e a k e k e )()()(1θ n b k e b k e ∂∂∂∂)()(0 Tn c k e c k e ⎥⎦⎤∂∂∂∂)()(1--------++-+∂∂=∂∂)()1()()()1()([)(101n k u b k u b k u b n k y a k y a k y a a k e n n i i )]()1(1n k e c k e c n ----in i i a n k e c a k e c a k e c i k y ∂-∂--∂-∂-∂-∂--=)()2()1()(21 （2.28） 即i nj j i a j k e c i k y a k e ∂-∂--=∂∂∑=)()()(1（2.29） 同理可得i nj j i b j k e c i k u b k e ∂-∂---=∂∂∑=)()()(1 （2.30） in j j i c j k e c i k e c k e ∂-∂---=∂∂∑=)()()(1 （2.31） 将式（2.29）移项化简，有in j j i n j j i a j k e c a j k e c a k e i k y ∂-∂=∂-∂+∂∂=-∑∑==)()()()(01 （2.32）因为j z k e j k e -=-)()( （2.33）由)(j k e -求偏导，故iji a z k e a j k e ∂∂=∂-∂-)()( （2.34） 将（2.34）代入（2.32），所以j nj j i i j n j j i nj j z c a k e a z k e c a j k e c i k y -=-==∑∑∑∂∂=∂∂=∂-∂=-000)()()()( （2.35） n n z c z c z c ---+++= 1111)(所以得)()()(1i k y a k e z c i-=∂∂- （2.36） 同理可得（2.30）和（2.31）为)()()(1i k u b k e z c i--=∂∂- （2.37） )()()(1i k e c k e z c i--=∂∂- （2.38） 根据（2.36）构造公式)(])([)]([)(1i k y j j i k y a j i k e z c j-=---=∂--∂- （2.39）将其代入（2.36），可得ij a k e z c a j i k e z c ∂∂=∂--∂--)()()]([)(11 (2.40）消除)(1-z c 可得1)1()()(a i k e a j i k e a k e j i ∂+-∂=∂+-∂=∂∂ （2.41） 同理可得（2.37）和（2.38）式)()()(b i k e b j i k e b k e j i ∂-∂=∂+-∂=∂∂ （2.42）1)1()()(c i k e c j i k e c k e j i ∂+-∂=∂+-∂=∂∂ （2.43） 式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为0，可通过求解这些差分方程，分别求出e(k)关于n n c c b b a a ,,,,,10,1的全部偏导数，而这些偏导数分别为)}({k y ，)}({k u 和)}({k e 的线性函数。