用極大似然法進行參數(shù)估計

1、.專業(yè)整理 .北京工商大學系統(tǒng)辨識 課程上機實驗報告（2014 年秋季學期 ）專業(yè)名稱：控制工程上機題目：極大似然法進行參數(shù)估計專業(yè)班級：2015年1月. 學習幫手 .專業(yè)整理 .一 實驗目的通過實驗掌握極大似然法在系統(tǒng)參數(shù)辨識中的原理和應用。二 實驗原理1 極大似然原理設(shè)有離散隨機過程 Vk 與未知參數(shù)有關(guān) ，假定已知概率分布密度f (Vk) 。如果我們得到 n 個獨立的觀測值V1,V2 ,nf (V1)，f (V2), f (Vn)。,V ，則可得分布密度， 要求根據(jù)這些觀測值來估計未知參數(shù)，估計的準則是觀測值 V 的出現(xiàn)概率為最大。k為此 ，定義一個似然函數(shù)L(V1,V2, ,Vn)f

2、(V1 ) f (V2 )f (Vn )（1.1）上式的右邊是 n 個概率密度函數(shù)的連乘，似然函數(shù) L是的函數(shù) 。 如果 L 達到極大值 ， Vk 的出現(xiàn)概率為最大 。 因此 ，極大似然法的實質(zhì)就是求出使L 達到極大值的的估值。為了便于求，對式 （ 1.1）等號兩邊取對數(shù) ，則把連乘變成連加，即nln Lln f (Vi )（1.2）i 1由于對數(shù)函數(shù)是單調(diào)遞增函數(shù)，當 L 取極大值時 ， lnL 也同時取極大值 。 求式 （1.2）對的偏導數(shù) ，令偏導數(shù)為0，可得ln L0（ 1.3）解上式可得的極大似然估計ML 。2 系統(tǒng)參數(shù)的極大似然估計Newton-Raphson法實際上就是一種遞推算

3、法，可以用于在線辨識。不過它是一種依每 L 次觀測數(shù)據(jù)遞推一次的算法 ，現(xiàn)在我們討論的是每觀測一次數(shù)據(jù)就遞推計算一次參數(shù)估計值得算法 。 本質(zhì)上說 ，它只是一種近似的極大似然法 。設(shè)系統(tǒng)的差分方程為. 學習幫手 .專業(yè)整理 .(1)(k)b(z1)u(k)(k)（2.1）a zy式中a( z 1) 1a1 z 1.an z nb( z 1 )b0b1z 1.bn z n因為(k) 是相關(guān)隨機向量 ，故（ 2.1）可寫成(z1) ()(z1)(k)(z1 ) (k)（2.2）ay kbuc式中c( z 1)(k )( k)（2.3）c( z 1)1c1z 1cn z n（2.4）(k) 是均值為

4、0的高斯分布白噪聲序列。 多項式 a( z 1)， b( z 1) 和 c(z 1 ) 中的系數(shù)a, a,b ,b ,c ,c 和序列 (k ) 的均方差都是未知參數(shù) 。1,0n1n設(shè)待估參數(shù)aab0bnc1cnT（2.5）1n并設(shè) y(k ) 的預測值為y( k)a1 y( k 1)an y(k n) b0 u( k)bn u(k n)(1)()（2.6）c1 e kcn e kn式中 e(ki) 為預測誤差 ； ai， bi， ci 為 a ， b ， c 的估值 。 預測誤差可表示為iiinne( k)y(k )y( k)y(k)aiy(ki)biu( ki )i 1i 0n1n1nci

5、 e(ki )1znz) y(k)(b01zbnz)u(k)(1 aabi1(c1z 1c2z 2cn zn )e( k)（2.7）或者(1 c1 z 1cn z n )e( k) = (1 a1 z 1an z n ) y(k )(b0b1 z 1bn z n )u(k )（2.8）因此預測誤差e(k )滿足關(guān)系式(z1 ) ()(z1)(k)(z1 )u(k)（2.9）ce kayb式中a( z 1) 1 a1 z 1an zb( z 1) b0 b1 z 1bn zc( z 1) 1 c1 z 1cn znnn假定預測誤差 e(k) 服從均值為 0 的高斯分布 ，并設(shè)序列e( k) 具有

6、相同的方差2。因. 學習幫手 .專業(yè)整理 .為 e(k) 與 c(z1 ) ， a( z1) 和 b( z 1) 有關(guān) ，所以2 是被估參數(shù)把式 （ 2.9）寫成(z1) ()(z1)(k)(z1 )(k)ce kaybue(k)y(k )a1 y(k1)an y(kn) b0 u(k 1)bnu(k n)c1e(k 1)cn (k n), k n 1,n或?qū)懗蒼nn的函數(shù) 。 為了書寫方便，（ 2.10 ）b1u(k1)2,（2.11 ）e(k )y( k)ai y( k i)bi u(k i)cie(k i )（2.12 ）i 1i 0i 1令 k=n+1,n+2,n+N, 可得 e(k)

7、 的 N 個方程式 ，把這 N 個方程式寫成向量-矩陣形式eNYN N（ 2.13 ）式中a1y(n1)e(n1)y(n2)e(n2)anYN, eN,b0y(nN )e(nN )bny( n)y(1)u( n1)u(1)e( n)e(1)y( n1)y( 2)u( n2)u(2)e( n1)e(2)Ny( nN1)y( N )u( n N )u(N )e( nN 1)e(N )因為已假定e(k ) 是均值為0 的高斯噪聲序列 ，高斯噪聲序列的概率密度函數(shù)為f11 exp1 2 ( y m)2 (22 ) 22（ 2.14 ）式中 y 為觀測值 ，2 和 m 為 y 的方差和均值 ，那么f11

8、exp1 2e2 (k)（2.15 ）(22 )22對于 e(k) 符合高斯噪聲序列的極大似然函數(shù)為L(YN ,)L e(n1),e( n2), e(nN ) f e(n1) f e(n2) f e(nN ) 1Nexp1 e2 (n1)e2 ( n2)e2 ( nN )1Nexp(1eNT eN )2) 222( 22222(2)（2.16 ）. 學習幫手 .專業(yè)整理 .或L(YN , )1ex p (YN)T (YN)（2.17 ）N22(22 ) 2對上式 （ 2.17 ）等號兩邊取對數(shù)得ln L (YN , ) ln1Nln exp(1eNTeN )N ln 2N ln221eNT e

9、N( 22) 222222（2.18 ）或?qū)憺镹 ln 2N ln1n Nln L (YN , )2e2 ( k)222 2k n 1求 ln L(Y, ) 對2 的偏導數(shù) ，令其等于0，可得Nln L (YN ,)N1nN2( k)02224e2k n 1則21 n N e2 (k)2 1 n N e2 (k)2 JN k n 1N 2 k n 1N式中（ 2.19 ）（ 2.20 ）（ 2.21 ）J1 nN2(k )（ 2.22 ）2 ken 12 越小越好 ，因為當方差2 最小時 ， e2 ( k) 最小，即殘差最小 。因此希望2 的估值取最小22 min J（ 2.23 ）N因為式

10、（ 2.10 ）可理解為預測模型，而 e(k) 可看做預測誤差 。因此使式 （ 2.22 ）最小就是使誤差的平方之和最小，即使對概率密度不作任何假設(shè)，這樣的準則也是有意義的。因此可按 J 最小來求 a1, , a,b0 ,bn ,c1,cn 的估計值 。由于 e(k)式參數(shù) a, a, b, b , c ,c 的線性函數(shù) ，因此 J 是這些參數(shù)的二次型函1,0n1n數(shù)。求使 ln L(YN, )最大的，等價于在式 （ 2.10 ）的約束條件下求使 J 為最小 。由于J 對 ci 是非線性的，因而求 J 的極小值問題并不好解 ，只能用迭代方法求解。求 J 極小值的常用迭代算法有拉格朗日乘子法和牛

11、頓- 拉卜森法 。 下面介紹牛頓- 拉卜森法 。 整個迭代計算步驟如下 ：(1)確定初始的0 值。對于0 中的 a1, a, b0 ,bn 可按模型e(k )a(z 1 ) y(k)b( z 1)u( k)（ 2.24 ）. 學習幫手 .專業(yè)整理 .用最小二乘法來求，而對于0中的 c1,cn 可先假定一些值 。(2)計算預測誤差e(k)y(k)y( k)（ 2.25）給出J1 n Ne2 (k)2 kn1并計算21n N2(k )N ken1J2J2 ,有(3)計算 J的梯度和海賽矩陣JnNe(k)e(k )kn 1式中T（ 2.26 ）（ 2.27 ）e(k )e(k)e(k)e( k)e(

12、k)e(k)e(k)a1anb0bnc1cne(k) y(k ) a1y( k 1)an y(k n) b0u(k ) b1u(k 1)bnu(k n)aiaic1e(k 1)cne(kn)y(k i ) c1e( k 1)e(k 2)e(k n)（ 2.28 ）c2cnaiaiai即e(k)ne( k j )c jy( k i )（2.29）aij 1ai同理可得. 學習幫手 .專業(yè)整理 .e(k)ne(kj )i )c ju(kbibij1e(k)ne(kj )i)c je(kcicij1將式 （ 2.29 ）移項化簡 ，有e(k )ne( kj)ne(k j )cjy(k i )aicj

13、aij 1j 0ai因為e(kj )e(k) zj由 e(kj ) 求偏導 ，故e(kj )e(k) z jaiai將（ 2.34 ）代入 （ 2.32 ）， 所以ne(k j )ne( k) z je( k)nc j zjy(k i)cjc jaij 0aij 0aij 0c( z 1) 1 c1 z 1cn z n所以得c( z 1) e(k)y(k i )ai同理可得 （ 2.30 ）和（ 2.31 ）為c( z 1)e(k)u(k i )bic( z1)e( k)(i)cie k根據(jù) （ 2.36 ）構(gòu)造公式c( z 1) e k (ij )yk (ij ) j y(k i )a j（

14、 2.30 ）（ 2.31 ）（ 2.32 ）（ 2.33 ）（ 2.34 ）（ 2.35 ）（ 2.36 ）（ 2.37 ）（ 2.38 ）（ 2.39 ）. 學習幫手 .專業(yè)整理 .將其代入 （ 2.36 ），可得c( z 1 ) e k(ij )c( z 1) e(k)(2.40 ）a jai消除 c( z 1) 可得e(k)e(k i j )e(k i 1)）aiaj（ 2.41a1同理可得 （ 2.37 ）和（ 2.38 ）式e(k)e( kij )e(ki )（2.42 ）bib jb0e(k)e(kij )e(ki 1)（ 2.43 ）cic jc1式（ 2.29 ）、 式（ 2

15、.30 ）和式 （ 2.31 ）均為差分方程 ，這些差分方程的初始條件為0，可通過求解這些差分方程，分別求出e(k) 關(guān)于 a1, ,a,b0 ,bn ,c1,cn 的全部偏導數(shù) ，而這些偏導數(shù)分別為 y(k ) ， u(k ) 和 e(k) 的線性函數(shù) 。 下面求關(guān)于的二階偏導數(shù) ，即2T2n Ne(k)e(k )nNJe(k )e(k)（2.44 ）22kn 1kn 1當接近于真值時， e(k) 接近于0。 在這種情況下 ，式（ 2.44 ）等號右邊第 2項接近于 0，2 J 可近似表示為22Jn Ne( k)Te(k )（2.45 ）2k n 1則利用式 （ 2.45 ）計算2 J比較簡

16、單 。2(4)按牛頓 - 拉卜森計算的新估值1 ，有10(2J)1J（2.46 ）20. 學習幫手 .專業(yè)整理 .重復 (2)至 (4) 的計算步驟 ，經(jīng)過 r 次迭代計算之后可得r ，近一步迭代計算可得r 1r(2J)1J（ 2.47 ）2r如果22r 1r4102r（2.48 ）則可停止計算 ，否則繼續(xù)迭代計算 。式（ 2.48 ）表明 ，當殘差方差的計算誤差小于0.01% 時就停止計算 。 這一方法即使在噪聲比較大的情況也能得到較好的估計值。三 實驗內(nèi)容設(shè) SISO 系統(tǒng)差分方程為y( k)a1 y(k1)a2 y(k2)b1u(k1)b2u(k2)(k )其中極大似然估計法默認e(k)

17、 為：e(k )C (z 1 )(k)(k )c1( k1) c2 ( k 2)辨識參數(shù)向量為 a1a2b1b2c1c2 T式中 ， ( k ) 為噪聲方差各異的白噪聲或有色噪聲；u(k) 和 y(k) 表示系統(tǒng)的輸入輸出變量 。. 學習幫手 .專業(yè)整理 .四 實驗流程圖五 代碼實現(xiàn)程序如下 ：U=1.147 0.201 -0.787 -1.584 -1.052 0.866 1.152 1.573 0.626.0.433-0.9580.810 -0.044 0.947 -1.474-0.719 -0.086 1.099.1.4501.151 0.485 1.633 0.043 1.326 1.

18、706 -0.340 0.890.0.433-1.177-0.390-0.9821.435-0.119 -0.769 -0.899 0.882.-1.008 -0.844 0.628-0.6791.5411.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.2031.204 1.691 -1.235-1.228 -1.267 0.309. 學習幫手 .專業(yè)整理 .0.0430.043 1.4611.585 0.552 -0.601 -0.3190.744

19、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.2880.8461.3250.7230.7130.643 0.463.0.7861.161 0.850-1.349 -0.5961.5120.795-0.713 0.453.-1.604 0.889-0.9380.0560.829-0.981 -1.232 1.327 -0.681.0.114-1.1351.28

20、4 -1.2010.7580.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.8030.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.9870.167-1.445 0.630 1.255 0.311 -1.726.0.9751.718 1.3601.667 1.111 1.018 0.078

21、 -1.665 -0.760.1.184-0.6140.994 -0.0890.9471.706-0.395 1.222 -1.351.0.2311.425 0.114-0.689 -0.7041.0700.2621.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.6420.805-2.088 0.946-0.0391.984-2.545-1.727 -

22、0.231.2.4403.5832.9151.443 3.598 0.702 2.638 3.611 -0.168.1.7320.6662.377-0.554-2.0882.6980.189 -1.633 -2.010.1.716-1.641 -1.885 1.061 -0.968 2.911 3.088-1.629 -1.533. 學習幫手 .專業(yè)整理 .3.0300.614-1.483-1.029-1.948-1.066-0.113 -2.144 -2.626.0.134-3.043-1.3410.3382.7023.813 -1.924-2.813-1.795.3.0021.0271.0

23、27 2.755 3.584 1.737 -0.837 -0.617 1.703.2.045-2.886-0.542-2.991-1.8593.0450.068-0.375 0.451.1.0360.153-0.4742.512 -2.681-0.954-0.3070.628-0.270.-0.277 0.983-0.2880.8461.3250.723 1.750 1.401 1.340.0.9161.3962.446 2.103 2.432 -1.486 3.031 2.373 -0.763.-0.752 -3.207 1.385-1.642-0.1181.756-1.613 -1.690

24、 2.136.-1.136 -0.005 -2.210 2.331-2.2040.9831.347-1.691 0.595.1.809-2.204-2.330-0.4541.2902.080-1.990-0.770 1.240.-0.252 3.137-2.3791.2061.221-1.9772.471-1.6801.148.1.8160.055-1.8560.269 -1.323-2.4861.9580.823 2.481.2.2093.167-0.762-2.225-0.123-2.7861.0262.8431.071.-3.317 1.5143.8073.388 3.683 -1.93

25、5 -1.9350.309 -3.390.-2.124 2.192-0.855-1.6560.0161.8043.774-0.0592.371.-2.322 -0.032 2.6320.565-1.460-1.8391.9170.8653.180.3.261-2.755-0.536-1.171-0.905-3.303 -0.834 2.490 3.039.0.1341.901% 輸出數(shù)據(jù)na=2;nb=1;nc=2;d=1;nn=max(na,nc);L=size(Y,2);. 學習幫手 .專業(yè)整理 .xiek=zeros(nc,1); %白噪聲估計初值yfk=zeros(nn,1); %yf

26、(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=(

27、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ù). 學習幫手 .專業(yè)整理 .thetae_1=thetae(:,k);for i=nc:-1:2xiek(i)=xiek(i-1);xiefk(i)=xiefk(i-1);endxiek(1)=xie;xie

28、fk(1)=xief;for i=nn:-1:2yfk(i)=yfk(i-1);ufk(i)=ufk(i-1);endyfk(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,:);. 學習幫手 .專業(yè)整理 .xlabel(k); ylabel(參數(shù)估計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(參數(shù)估計c);legend(c_1,c_