用极大似然法进行参数估计.docx
《用极大似然法进行参数估计.docx》由会员分享，可在线阅读，更多相关《用极大似然法进行参数估计.docx（65页珍藏版）》请在冰豆网上搜索。
用极大似然法进行参数估计
.专业整理.
北京工商大学
《系统辨识》课程
上机实验报告
（2014年秋季学期）
专业名称：
控制工程
上机题目：
极大似然法进行参数估计
专业班级：
2015年1月
.学习帮手.
.专业整理.
一实验目的
通过实验掌握极大似然法在系统参数辨识中的原理和应用。
二实验原理
1极大似然原理
设有离散随机过程{Vk}与未知参数
有关，假定已知概率分布密度
f（Vk
）。
如果我
们得到n个独立的观测值
V1,V2,
n
f（V1
）
，
f（V2
）
f（Vn
）
。
V，则可得分布密度
，
要求根据这些观测值来估计未知参数
，估计的准则是观测值
{{V
}}的出现概率为最大
。
k
为此，定义一个似然函数
L（V1,V2,,Vn
）
f（V1）f（V2）
f（Vn）
（1.1）
上式的右边是n个概率密度函数的连乘
，似然函数L是
的函数。
如果L达到极大值，{Vk}
的出现概率为最大。
因此，极大似然法的实质就是求出使
L达到极大值的
的估值
。
为
了便于求
，对式（1.1）等号两边取对数，则把连乘变成连加
，即
n
lnL
lnf（Vi）
（1.2）
i1
由于对数函数是单调递增函数
，当L取极大值时，lnL也同时取极大值。
求式（1.2）对
的偏导数，令偏导数为
0，可得
lnL
0
（1.3）
解上式可得的极大似然估计ML。
2系统参数的极大似然估计
Newton-Raphson法实际上就是一种递推算法，可以用于在线辨识。
不过它是一种依
每L次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。
本质上说，它只是一种近似的极大似然法。
设系统的差分方程为
.学习帮手.
.专业整理.
（
1）
（
k
）
b
（
z
1）
u
（
k
）
（
k
）
（2.1）
az
y
式中
a（z1）1
a1z1
...
anzn
b（z1）
b0
b1z1
...
bnzn
因为
（k）是相关随机向量，故（2.1）可写成
（
z
1）（
）
（
z
1）
（
k
）
（
z
1）（
k
）
（
2.2
）
a
yk
b
u
c
式中
c（z1）
（k）
（k）
（2.3）
c（z1）
1
c1z1
cnzn
（2.4）
（k）是均值为
0
的高斯分布白噪声序列
。
多项式a（z1）
，b（z1）和c（z1）中的系数
a
a,b,
b,c,
c和序列{（k）}的均方差
都是未知参数。
1,
0
n
1
n
设待估参数
[a
a
b0
bn
c1
cn
T
（2.5）
1
n
并设y（k）的预测值为
y（k）
a1y（k1）
any（kn）b0u（k）
bnu（kn）
（
1）
（
）
（
2.6
）
c1ek
cnek
n
式中e（k
i）为预测误差；ai
，bi
，ci为a，b，c的估值。
预测误差可表示为
i
i
i
n
n
e（k）
y（k）
y（k）
y（k）
ai
y（k
i）
biu（k
i）
i1
i0
n
1
n
1
n
cie（k
i）
1
z
n
z
）y（k）
（b
0
1
z
b
n
z
）u（k）
（1a
a
b
i
1
（c1
z1
c2
z2
cnz
n）e（k）
（2.7）
或者
（1c1z1
cnzn）e（k）=（1a1z1
anzn）y（k）
（b0
b1z1
bnzn）u（k）
（2.8）
因此预测误差
e（k）
满足关系式
（
z
1）（
）
（
z
1）
（
k
）
（
z
1）
u
（
k
）
（2.9）
c
ek
a
y
b
式中
a（z1）1a1z1
anz
b（z1）b0b1z1
bnz
c（z1）1c1z1
cnz
n
n
n
假定预测误差e（k）服从均值为0的高斯分布，并设序列
e（k）具有相同的方差
2
。
因
.学习帮手.
.专业整理.
为e（k）与c（z
1），a（z
1）和b（z1）有关，所以
2是被估参数
把式（2.9）写成
（
z
1）（
）
（
z
1）
（
k
）
（
z
1）
（
k
）
c
ek
a
y
b
u
e（k）
y（k）
a1y（k
1）
any（k
n）b0u（k1）
bnu（kn）
c1e（k1）
cn（kn）,kn1,n
或写成
n
n
n
的函数。
为了书写方便，
（2.10）
b1u（k1）
2,（2.11）
e（k）
y（k）
aiy（ki）biu（ki）
cie（ki）
（2.12）
i1
i0
i1
令k=n+1,n+2,
n+N,可得e（k）的N个方程式，把这N个方程式写成向量
-矩阵形式
eNYNN
（2.13）
式中
a1
y（n
1）
e（n
1）
y（n
2）
e（n
2）
an
YN
eN
b0
y（n
N）
e（n
N）
bn
y（n）
y
（1）
u（n
1）
u
（1）
e（n）
e
（1）
y（n
1）
y
（2）
u（n
2）
u
（2）
e（n
1）
e
（2）
N
y（n
N
1）
y（N）
u（nN）
u（N）
e（n
N1）
e（N）
因为已假定
e（k）是均值为
0的高斯噪声序列，高斯噪声序列的概率密度函数为
f
1
1exp[
12（ym）2]
（2
2）2
2
（2.14）
式中y为观测值，
2和m为y的方差和均值，那么
f
1
1
exp[
12
e2（k）]
（2.15）
（2
2）2
2
对于e（k）符合高斯噪声序列的极大似然函数为
L（YN,
）
L[e（n
1）,e（n
2）,
e（n
N）]
f[e（n
1）]f[e（n
2）
]f[e（n
N）]
1
N
exp{
1
[e2（n
1）
e2（n
2）
e2（n
N）]}
1
N
exp（
1
eNTeN）
2
）2
2
2
（2
2
2
2
2
（2
）
（2.16）
.学习帮手.
.专业整理.
或
L（YN,）
1
exp[
（YN
）T（YN
）
]
（2.17）
N
2
2
（2
2）2
对上式（2.17）等号两边取对数得
lnL（YN,）ln
1
N
lnexp（
1
eNTeN）
Nln2
Nln
2
2
1
eNTeN
（2
2
）2
2
2
2
2
2
（2.18）
或写为
Nln2
Nln
1
nN
lnL（YN,）
2
e2（k）
2
2
22
kn1
求lnL（Y
）对
2的偏导数，令其等于
0，可得
N
lnL（YN,
）
N
1
n
N
2
（k）
0
2
2
2
4
e
2
kn1
则
2
1nNe2（k）
21nNe2（k）
2J
Nkn1
N2kn1
N
式中
（2.19）
（2.20）
（2.21）
J
1n
N
2
（k）
（2.22）
2k
e
n1
2越小越好，因为当方差
2最小时，e2（k）最小，即残差最小。
因此希望
2的估值取最
小
2
2minJ
（2.23）
N
因为式（2.10）可理解为预测模型
，而e（k）可看做预测误差。
因此使式（2.22）最小就
是使误差的平方之和最小
，即使对概率密度不作任何假设
，这样的准则也是有意义的
。
因
此可按J最小来求a1,,a,b0,
bn,c1,
cn的估计值。
由于e（k）式参数a
a,b
b,c,
c的线性函数，因此J是这些参数的二次型函
1,
0
n
1
n
数。
求使lnL（YN
）最大的
，等价于在式（2.10）的约束条件下求
使J为最小。
由于
J对ci是非线性的
，因而求J的极小值问题并不好解，只能用迭代方法求解
。
求J极小值的
常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。
下面介绍牛顿-拉卜森法。
整个迭代计
算步骤如下：
（1）确定初始的0值。
对于0中的a1,,a,b0,bn可按模型
e（k）a（z1）y（k）b（z1）u（k）（2.24）
.学习帮手.
.专业整理.
用最小二乘法来求
，而对于0
中的c1,
cn可先假定一些值。
（2）计算预测误差
e（k）
y（k）
y（k）
（2.25
）
给出
J
1nN
e2（k）
2k
n
1
并计算
2
1
nN
2
（k）
Nk
e
n
1
J
2J2,有
（3）计算J的梯度
和海赛矩阵
J
n
N
e（k）
e（k）
k
n1
式中
T
（2.26）
（2.27）
e（k）e（k）e（k）e（k）e（k）e（k）e（k）
a1anb0bnc1cn
e（k）
[y（k）a1y（k1）
any（kn）b0u（k）b1u（k1）
bnu（kn）
ai
ai
c1e（k1）
cne（k
n）]
y（ki）c1
e（k1）
e（k2）
e（kn）
（2.28）
c2
cn
ai
ai
ai
即
e（k）
n
e（kj）
cj
y（ki）
（2.29
）
ai
j1
ai
同理可得
.学习帮手.
.专业整理.
e（k）
n
e（k
j）
i）
cj
u（k
bi
bi
j
1
e（k）
n
e（k
j）
i）
cj
e（k
ci
ci
j
1
将式（2.29）移项化简，有
e（k）
n
e（k
j）
n
e（kj）
cj
y（ki）
ai
cj
ai
j1
j0
ai
因为
e（k
j）
e（k）z
j
由e（k
j）求偏导，故
e（k
j）
e（k）zj
ai
ai
将（2.34）代入（2.32），所以
n
e（kj）
n
e（k）zj
e（k）
n
cjz
j
y（ki）
cj
cj
ai
j0
ai
j0
ai
j0
c（z1）1c1z1
cnzn
所以得
c（z1）e（k）
y（ki）
ai
同理可得（2.30）和（2.31）为
c（z1）
e（k）
u（ki）
bi
c（z
1
）
e（k）
（
i
）
ci
ek
根据（2.36）构造公式
c（z1）e[k（i
j）]
y[k（i
j）j]y（ki）
aj
（2.30）
（2.31）
（2.32）
（2.33）
（2.34）
（2.35）
（2.36）
（2.37）
（2.38）
（2.39）
.学习帮手.
.专业整理.
将其代入（2.36），可得
c（z1）e[k（i
j）]
c（z1）e（k）
（2.40）
aj
ai
消除c（z1）可得
e（k）
e（kij）
e（ki1）
）
ai
aj
（2.41
a1
同理可得（2.37）和（2.38）式
e（k）
e（k
i
j）
e（k
i）
（2.42）
bi
bj
b0
e（k）
e（k
i
j）
e（k
i1）
（2.43）
ci
cj
c1
式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为
0，
可通过求解这些差分方程
，分别求出
e（k）关于a1,,a,b0,
bn,c1,
cn的全部偏导数，而这
些偏导数分别为{y（k）}，{u（k）}和{e（k）}的线性函数。
下面求关于
的二阶偏导数，即
2
T
2
nN
e（k）
e（k）
n
N
J
升级会员 