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

北京工商大学 《系统辨识》课程 上机实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 （2014年秋季学期） 专业名称 ：        控制工程                上机 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 目 ：      极大似然法进行参数估计    专业班级 ：                    2015年 1 月 一 实验目的 通过实验掌握极大似然法在系统参数辨识中的原理和应用。 二 实验原理 1 极大似然原理 设有离散随机过程 与未知参数 有关，假定已知概率分布密度 。如果我们得到n个独立的观测值 … ，则可得分布密度 ， ，…, 。要求根据这些观测值来估计未知参数 ，估计的准则是观测值{ }的出现概率为最大。为此，定义一个似然函数 （1.1）                      上式的右边是n个概率密度函数的连乘，似然函数L是 的函数。如果L达到极大值， 的出现概率为最大。因此，极大似然法的实质就是求出使L达到极大值的 的估值 。为了便于求 ，对式（1.1）等号两边取对数，则把连乘变成连加，即 （1.2） 由于对数函数是单调递增函数，当L取极大值时，lnL也同时取极大值。求式（1.2）对 的偏导数，令偏导数为0，可得 （1.3） 解上式可得 的极大似然估计 。 2 系统参数的极大似然估计 Newton-Raphson法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。本质上说，它只是一种近似的极大似然法。 设系统的差分方程为 （2.1） 式中 因为 是相关随机向量，故（2.1）可写成 （2.2） 式中 （2.3） （2.4） 是均值为0的高斯分布白噪声序列。多项式 ， 和 中的系数 和序列 的均方差 都是未知参数。 设待估参数 （2.5） 并设 的预测值为 （2.6） 式中 为预测误差； ， ， 为 ， ， 的估值。预测误差可表示为 （2.7） 或者 = （2.8） 因此预测误差 满足关系式 （2.9） 式中 假定预测误差 服从均值为0的高斯分布，并设序列 具有相同的方差 。因为 与 ， 和 有关，所以 是被估参数 的函数。为了 书 关于书的成语关于读书的排比句社区图书漂流公约怎么写关于读书的小报汉书pdf 写方便，把式（2.9）写成 （2.10） （2.11） 或写成 （2.12） 令k=n+1,n+2,…,n+N,可得 的N个方程式，把这N个方程式写成向量-矩阵形式 （2.13） 式中 , , 因为已假定 是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为 （2.14） 式中y为观测值， 和m为y的方差和均值，那么 （2.15） 对于 符合高斯噪声序列的极大似然函数为 （2.16） 或 （2.17） 对上式（2.17）等号两边取对数得 （2.18）        或写为 （2.19） 求 对 的偏导数，令其等于0，可得 （2.20） 则 （2.21） 式中 （2.22） 越小越好，因为当方差 最小时， 最小，即残差最小。因此希望 的估值取最小 （2.23） 因为式（2.10）可理解为预测模型，而e(k)可看做预测误差。因此使式（2.22）最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因此可按J最小来求 的估计值。 由于e(k)式参数 的线性函数，因此J是这些参数的二次型函数。求使 最大的 ，等价于在式（2.10）的约束条件下求 使J为最小。由于J对 是非线性的，因而求J的极小值问题并不好解，只能用迭代方法求解。求J极小值的常用迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下： (1)确定初始的 值。对于 中的 可按模型 （2.24） 用最小二乘法来求，而对于 中的 可先假定一些值。 (2)计算预测误差 （2.25） 给出  并计算                （2.26） (3)计算J的梯度 和海赛矩阵 ,有 （2.27） 式中  （2.28） 即 （2.29） 同理可得 （2.30） （2.31） 将式（2.29）移项化简，有 （2.32） 因为 （2.33） 由 求偏导，故 （2.34） 将（2.34）代入（2.32），所以 （2.35） 所以得 （2.36） 同理可得（2.30）和（2.31）为 （2.37） （2.38） 根据（2.36）构造公式 （2.39） 将其代入（2.36），可得 (2.40） 消除 可得 （2.41） 同理可得（2.37）和（2.38）式 （2.42） （2.43） 式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为0，可通过求解这些差分方程，分别求出e(k)关于 的全部偏导数，而这些偏导数分别为 ， 和 的线性函数。下面求关于 的二阶偏导数，即 （2.44）        当 接近于真值 时，e(k)接近于0。在这种情况下，式（2.44）等号右边第2项接近于0， 可近似表示为 （2.45） 则利用式（2.45）计算 比较简单。 (4)按牛顿-拉卜森计算 的新估值 ，有 （2.46） 重复(2)至(4)的计算步骤，经过r次迭代计算之后可得 ，近一步迭代计算可得 （2.47） 如果 （2.48） 则可停止计算，否则继续迭代计算。 式（2.48）表明，当残差方差的计算误差小于 时就停止计算。这一方法即使在噪声比较大的情况也能得到较好的估计值 。 三 实验内容 设SISO系统差分方程为 其中极大似然估计法默认 为： 辨识参数向量为 ＝         c1  c2] 式中， 为噪声方差各异的白噪声或有色噪声；u(k)和y(k)表示系统的输入输出变量。 四 实验 流程 快递问题件怎么处理流程河南自建厂房流程下载关于规范招聘需求审批流程制作流程表下载邮件下载流程设计 图 五 代码实现 程序如下： U=[1.147 0.201 -0.787 -1.584 -1.052 0.866 1.152 1.573 0.626... 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 1.099... 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890... 0.433 -1.177 -0.390 -0.982 1.435 -0.119 -0.769 -0.899 0.882... -1.008 -0.844 0.628 -0.679 1.541 1.375 -0.984 -0.582 1.609... 0.090 -0.813 -0.428 -0.848 -0.410 0.048 -1.099 -1.108 0.259... -1.627 -0.528 0.203 1.204 1.691 -1.235 -1.228 -1.267 0.309... 0.043 0.043 1.461 1.585 0.552 -0.601 -0.319 0.744 0.829... -1.626 -0.127 -1.578 -0.822 1.469 -0.379 -0.212 0.178 0.493... -0.056 -0.1294 1.228 -1.606 -0.382 -0.229 0.313 -0.161 -0.810... -0.277 0.983 -0.288 0.846 1.325 0.723 0.713 0.643 0.463... 0.786 1.161 0.850 -1.349 -0.596 1.512 0.795 -0.713 0.453... -1.604 0.889 -0.938 0.056 0.829 -0.981 -1.232 1.327 -0.681... 0.114 -1.135 1.284 -1.201 0.758 0.590 -1.007 0.390 0.836... -1.52 -1.053 -0.083 0.619 0.840 -1.258 -0.354 0.629 -0.242... 1.680 -1.236 -0.803 0.537 -1.100 1.417 -1.024 0.671 0.688... -0.123 -0.952 0.232 -0.793 -1.138 1.154 0.206 1.196 1.013... 1.518 -0.553 -0.987 0.167 -1.445 0.630 1.255 0.311 -1.726... 0.975 1.718 1.360 1.667 1.111 1.018 0.078 -1.665 -0.760... 1.184 -0.614 0.994 -0.089 0.947 1.706 -0.395 1.222 -1.351... 0.231 1.425 0.114 -0.689 -0.704 1.070 0.262 1.610 1.489... -1.602 0.020 -0.601 -0.020 -0.601 -0.235 1.245 1.226 -0.204... 0.926 -1.297] %输入数据 Y=[0.086 2.210 0.486 -1.812 -3.705 -2.688 1.577 2.883 3.705... 1.642 0.805 -2.088 0.946 -0.039 1.984 -2.545 -1.727 -0.231... 2.440 3.583 2.915 1.443 3.598 0.702 2.638 3.611 -0.168... 1.732 0.666 2.377 -0.554 -2.088 2.698 0.189 -1.633 -2.010... 1.716 -1.641 -1.885 1.061 -0.968 2.911 3.088 -1.629 -1.533... 3.030 0.614 -1.483 -1.029 -1.948 -1.066 -0.113 -2.144 -2.626... 0.134 -3.043 -1.341 0.338 2.702 3.813 -1.924 -2.813 -1.795... 3.002 1.027 1.027 2.755 3.584 1.737 -0.837 -0.617 1.703... 2.045 -2.886 -0.542 -2.991 -1.859 3.045 0.068 -0.375 0.451... 1.036 0.153 -0.474 2.512 -2.681 -0.954 -0.307 0.628 -0.270... -0.277 0.983 -0.288 0.846 1.325 0.723 1.750 1.401 1.340... 0.916 1.396 2.446 2.103 2.432 -1.486 3.031 2.373 -0.763... -0.752 -3.207 1.385 -1.642 -0.118 1.756 -1.613 -1.690 2.136... -1.136 -0.005 -2.210 2.331 -2.204 0.983 1.347 -1.691 0.595... 1.809 -2.204 -2.330 -0.454 1.290 2.080 -1.990 -0.770 1.240... -0.252 3.137 -2.379 1.206 1.221 -1.977 2.471 -1.680 1.148... 1.816 0.055 -1.856 0.269 -1.323 -2.486 1.958 0.823 2.481... 2.209 3.167 -0.762 -2.225 -0.123 -2.786 1.026 2.843 1.071... -3.317 1.514 3.807 3.388 3.683 -1.935 -1.935 0.309 -3.390... -2.124 2.192 -0.855 -1.656 0.016 1.804 3.774 -0.059 2.371... -2.322 -0.032 2.632 0.565 -1.460 -1.839 1.917 0.865 3.180... 3.261 -2.755 -0.536 -1.171 -0.905 -3.303 -0.834 2.490 3.039... 0.134 1.901]%输出数据 na=2;nb=1;nc=2;d=1; nn=max(na,nc); L=size(Y,2); xiek=zeros(nc,1); %白噪声估计初值 yfk=zeros(nn,1); %yf(k-i) ufk=zeros(nn,1); %uf(k-i) xiefk=zeros(nc,1); %vf(k-i) thetae_1=zeros(na+nb+1+nc,1); %参数估计初值 P=eye(na+nb+1+nc); for k=3:L %构造向量 phi=[-Y(k-1);-Y(k-2);U(k-1);U(k-2);xiek]; %组建h（k） xie=Y(k)-phi'*thetae_1; phif=[-yfk(1:na);ufk(d:d+nb);xiefk]; %递推极大似然参数估计算法 K=P*phif/(1+phif'*P*phif); thetae(:,k)=thetae_1+K*xie; P=(eye(na+nb+1+nc)-K*phif')*P;    yf=Y(k)-thetae(na+nb+2:na+nb+1+nc,k)'*yfk(1:nc); %yf(k) uf=U(k)-thetae(na+nb+2:na+nb+1+nc,k)'*ufk(1:nc); %uf(k) xief=xie-thetae(na+nb+2:na+nb+1+nc,k)'*xiefk(1:nc); %vf(k) %更新数据 thetae_1=thetae(:,k); for i=nc:-1:2 xiek(i)=xiek(i-1); xiefk(i)=xiefk(i-1); end xiek(1)=xie; xiefk(1)=xief; for i=nn:-1:2 yfk(i)=yfk(i-1); ufk(i)=ufk(i-1); end yfk(1)=yf; ufk(1)=uf; end thetae_1 figure(1) plot([1:L],thetae(1:na,:)); xlabel('k'); ylabel('参数估计a'); legend('a_1','a_2'); axis([0 L -2 2]); figure(2) plot([1:L],thetae(na+1:na+nb+1,:)); xlabel('k'); ylabel('参数估计b'); legend('b_1','b_2'); axis([0 L -1.5 2]); figure(3) plot([1:L],thetae(na+nb+2:na+nb+nc+1,:)); xlabel('k'); ylabel('参数估计c'); legend('c_1','c_2'); axis([0 L -2 2]); 六 实验结果 实验结果如下所示。 thetae_1 = -0.1222 0.0365 1.7962 0.0004 0.0655 -0.0018 可知，估计值a1=-0.1222，a2=0.0365，b1=1.7962，b2= 0.0004，c1= 0.0655，c2= -0.0018 图1 A参数波形变化情况 图2 B参数波形变化情况 图3 C参数波形变化情况 图10-11 图10-12有162个数据，结果如下： thetae_1 = -0.7095 0.0489 1.8688 -0.9850 -0.1768 -0.0005 可知，估计值a1=-0.7095，a2=0.0489，b1= 1.8688，b2= -0.9850，c1= -0.1768，c2= -0.0005 图1 A参数波形变化情况 图2 B参数波形变化情况 图3 C参数波形变化情况 六 实验 总结 初级经济法重点总结下载党员个人总结TXt高中句型全总结.doc高中句型全总结.doc理论力学知识点总结pdf 通过此次实验，对极大似然法的原理和方法进行了研究和对其matlab的仿真，可发现其离线辨识效果很好，但是其在线辨识的计算量很大，实验数据不同，出现的效果也不一样，递推的极大似然法呈现病态，所得参数估计取值很不稳定，以至于不能采用。