用极大似然法进行参数估计

北京工商大学

《系统辨识》课程上机实验报告（2014年秋季学期）

专业名称：控制工程

上机题目：极大似然法进行参数估计专业班级：

2015年1 月

一 实验目的

通过实验掌握极大似然法在系统参数辨识中的原理和应用。

二 实验原理

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）等号两边取对数，则把连乘变成连加，即 ∑==

n

i i

V f L 1

)(ln ln θ （1.2）

由于对数函数是单调递增函数，当L 取极大值时，lnL 也同时取极大值。求式（1.2）对

θ的偏导数，令偏导数为0，可得

0ln =∂∂θL

（1.3）

解上式可得θ的极大似然估计ML ∧

θ。

2 系统参数的极大似然估计

Newton-Raphson 法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L 次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。本质上说，它只是一种近似的极大似然法。

设系统的差分方程为

)()()()()(1

1

k k u z b k y z a ξ+=-- （2.1）

式中

111()1...n

n a z a z a z ---=+++

1101()...n

n b z b b z b z

---=+++

因为)(k ξ是相关随机向量，故（2.1）可写成

)()()()()()(1

1

1

k z c k u z b k y z a ε---+= （2.2）

式中

)()()(1

k k z c ξε=- （2.3）

n

n 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 ]T

n 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 的估值。预测误差可表示为

+-+-⎢⎣⎡--=-=∑∑=∧

=∧

∧)()()()()()(01

i k u b i k y a k y k y k y k e n i i n i i

-+++-+++=⎥⎦

⎤--∧-∧∧-∧-∧=∧∑)()()()1()(1

10111k u z b z b b k y z a z a i k e c n n n n n i i )()(2

21

1k e z c z

c z c n n -∧

-∧

-∧+++ （2.7）

或者

)()1(11k e z c z c n

n -∧

-∧

+++ =-+++-∧

-∧

)()1(1

1k y z a z a n

n

)()(1

10k u z b z b b n

n -∧

-∧∧+++ （2.8）

因此预测误差{})(k e 满足关系式

)()()()()()(1

1

1

k u z b k y z a k e z c -∧

-∧

-∧

-= （2.9）

式中

n n z a z a z a -∧

-∧

-∧

+++= 1

11

1)( n n z b z b b z b -∧

-∧

∧-∧

+++= 1

101

)( n n z c z c z c -∧

-∧

-∧

+++= 1

11

1)(

假定预测误差)(k e 服从均值为0的高斯分布，并设序列{})(k e 具有相同的方差2

σ。因

为{})(k e 与)(1

-∧

z c ，)(1

-∧

z a 和)(1

-∧

z b 有关，所以2

σ是被估参数θ的函数。为了书写方便，

把式（2.9）写成

)()()()()()(1

1

1

k 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） 或写成

)()()()()(1

1

i k e c i k u b i k y a k y k e n

i i

n

i i

n

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(122

2

12m y f --

=

σ

πσ （2.14）

式中y 为观测值，2

σ和m 为y 的方差和均值，那么

)](21ex p[)

2(122

2

12k e f σ

πσ-

=

（2.15）

对于)(k e 符合高斯噪声序列的极大似然函数为

)

21exp()

2(1)]}()2()1([21exp{)

2(1])([])2([])1([])(,),2(),1([),(2

2

22222

2

2N T N N

N N 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）