用極大似然法進(jìn)行參數(shù)估計.doc

.北京工商大學(xué)系統(tǒng)建模與辨識課程上機實驗報告（2016年秋季學(xué)期）專業(yè)名稱 ： 控制工程 上機題目 ： 用極大似然法進(jìn)行參數(shù)估計 專業(yè)班級 ： 計研3班 學(xué)生姓名 ： 王 瑤 吳 超 學(xué) 號 ： 10011316259 10011316260 指導(dǎo)教師 ： 劉翠玲 2017 年 1 月一 實驗?zāi)康耐ㄟ^實驗掌握極大似然法在系統(tǒng)參數(shù)辨識中的原理和應(yīng)用。二 實驗原理1 極大似然原理設(shè)有離散隨機過程與未知參數(shù)有關(guān)，假定已知概率分布密度。如果我們得到n個獨立的觀測值，則可得分布密度，,。要求根據(jù)這些觀測值來估計未知參數(shù)，估計的準(zhǔn)則是觀測值的出現(xiàn)概率為最大。為此，定義一個似然函數(shù) （1.1） 上式的右邊是n個概率密度函數(shù)的連乘，似然函數(shù)L是的函數(shù)。如果L達(dá)到極大值，的出現(xiàn)概率為最大。因此，極大似然法的實質(zhì)就是求出使L達(dá)到極大值的的估值。為了便于求，對式（1.1）等號兩邊取對數(shù)，則把連乘變成連加，即 （1.2）由于對數(shù)函數(shù)是單調(diào)遞增函數(shù)，當(dāng)L取極大值時，lnL也同時取極大值。求式（1.2）對的偏導(dǎo)數(shù)，令偏導(dǎo)數(shù)為0，可得 （1.3）解上式可得的極大似然估計。 2 系統(tǒng)參數(shù)的極大似然估計Newton-Raphson法實際上就是一種遞推算法，可以用于在線辨識。不過它是一種依每L次觀測數(shù)據(jù)遞推一次的算法，現(xiàn)在我們討論的是每觀測一次數(shù)據(jù)就遞推計算一次參數(shù)估計值得算法。本質(zhì)上說，它只是一種近似的極大似然法。設(shè)系統(tǒng)的差分方程為 （2.1）式中 因為是相關(guān)隨機向量，故（2.1）可寫成 （2.2）式中 （2.3） （2.4）是均值為0的高斯分布白噪聲序列。多項式，和中的系數(shù)和序列的均方差都是未知參數(shù)。設(shè)待估參數(shù) （2.5）并設(shè)的預(yù)測值為 （2.6）式中為預(yù)測誤差；，為，的估值。預(yù)測誤差可表示為 （2.7）或者 = （2.8）因此預(yù)測誤差滿足關(guān)系式 （2.9）式中假定預(yù)測誤差服從均值為0的高斯分布，并設(shè)序列具有相同的方差。因為與，和有關(guān)，所以是被估參數(shù)的函數(shù)。為了書寫方便，把式（2.9）寫成 （2.10） （2.11）或?qū)懗?（2.12）令k=n+1,n+2,n+N,可得的N個方程式，把這N個方程式寫成向量-矩陣形式 （2.13）式中 , , 因為已假定是均值為0的高斯噪聲序列，高斯噪聲序列的概率密度函數(shù)為 （2.14）式中y為觀測值，和m為y的方差和均值，那么 （2.15）對于符合高斯噪聲序列的極大似然函數(shù)為 （2.16）或 （2.17）對上式（2.17）等號兩邊取對數(shù)得 （2.18） 或?qū)憺?（2.19）求對的偏導(dǎo)數(shù)，令其等于0，可得 （2.20）則 （2.21）式中 （2.22）越小越好，因為當(dāng)方差最小時，最小，即殘差最小。因此希望的估值取最小 （2.23）因為式（2.10）可理解為預(yù)測模型，而e(k)可看做預(yù)測誤差。因此使式（2.22）最小就是使誤差的平方之和最小，即使對概率密度不作任何假設(shè)，這樣的準(zhǔn)則也是有意義的。因此可按J最小來求的估計值。由于e(k)式參數(shù)的線性函數(shù)，因此J是這些參數(shù)的二次型函數(shù)。求使最大的，等價于在式（2.10）的約束條件下求使J為最小。由于J對是非線性的，因而求J的極小值問題并不好解，只能用迭代方法求解。求J極小值的常用迭代算法有拉格朗日乘子法和牛頓-拉卜森法。下面介紹牛頓-拉卜森法。整個迭代計算步驟如下：(1)確定初始的值。對于中的可按模型 （2.24）用最小二乘法來求，而對于中的可先假定一些值。(2)計算預(yù)測誤差 （2.25）給出 并計算 （2.26）(3)計算J的梯度 和海賽矩陣 ,有 （2.27）式中 （2.28）即 （2.29）同理可得 （2.30） （2.31）將式（2.29）移項化簡，有 （2.32）因為 （2.33）由求偏導(dǎo)，故 （2.34）將（2.34）代入（2.32），所以 （2.35）所以得 （2.36）同理可得（2.30）和（2.31）為 （2.37） （2.38）根據(jù)（2.36）構(gòu)造公式 （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)關(guān)于的全部偏導(dǎo)數(shù)，而這些偏導(dǎo)數(shù)分別為，和的線性函數(shù)。下面求關(guān)于的二階偏導(dǎo)數(shù)，即 （2.44） 當(dāng)接近于真值時，e(k)接近于0。在這種情況下，式（2.44）等號右邊第2項接近于0，可近似表示為 （2.45）則利用式（2.45）計算比較簡單。(4)按牛頓-拉卜森計算的新估值，有 （2.46）重復(fù)(2)至(4)的計算步驟，經(jīng)過r次迭代計算之后可得，近一步迭代計算可得 （2.47）如果 （2.48）則可停止計算，否則繼續(xù)迭代計算。式（2.48）表明，當(dāng)殘差方差的計算誤差小于時就停止計算。這一方法即使在噪聲比較大的情況也能得到較好的估計值。三 實驗內(nèi)容設(shè)SISO系統(tǒng)差分方程為 其中極大似然估計法默認(rèn)為： 辨識參數(shù)向量為 c1 c2式中，為噪聲方差各異的白噪聲或有色噪聲；u(k)和y(k)表示系統(tǒng)的輸入輸出變量。四 實驗流程圖五 代碼實現(xiàn)程序如下：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 %輸入數(shù)據(jù)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%輸出數(shù)據(jù)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); %參數(shù)估計初值P=eye(na+nb+1+nc);for k=3:L %構(gòu)造向量 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; %遞推極大似然參數(shù)估計算法 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) %更新數(shù)據(jù) 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;endthetae_1figure(1)plot(1:L,thetae(1:na,:);xlabel(k); ylabel(參數(shù)估計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(