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

1、北京系統(tǒng)辨識課程二機(jī)實(shí)驗(yàn)報(bào)告控制工程（2014年秋季學(xué)期）專業(yè)名稱：上機(jī)題目：極大似然法進(jìn)行參數(shù)估計(jì)專業(yè)班級：2015年1月實(shí)驗(yàn)?zāi)康耐ㄟ^實(shí)驗(yàn)掌握極大似然法在系統(tǒng)參數(shù)辨識中的原理和應(yīng)用。實(shí)驗(yàn)原理1極大似然原理設(shè)有離散隨機(jī)過程Vk與未知參數(shù) 有關(guān)，假定已知概率分布密度 fM )。如果我們Vk的出現(xiàn)概率為最大。為此,得到n個(gè)獨(dú)立的觀測值V1,V2,Vn，則可得分布密度f(V|) ,f(V2)，,f(Vn)。要求根據(jù)這些觀測值來估計(jì)未知參數(shù)，估計(jì)的準(zhǔn)則是觀測值定義一個(gè)似然函數(shù)L(Vi,V2, ,Vn ) f(Vi )f(V2 )f(Vn )上式的右邊是n個(gè)概率密度函數(shù)的連乘，似然函數(shù)L是(1.1)的函

2、數(shù)。如果L達(dá)到極大值，VkL達(dá)到極大值的 的估值。為了的出現(xiàn)概率為最大。因此，極大似然法的實(shí)質(zhì)就是求出使便于求，對式(1.1 )等號兩邊取對數(shù)，則把連乘變成連加，即In L由于對數(shù)函數(shù)是單調(diào)遞增函數(shù)，當(dāng) 對 的偏導(dǎo)數(shù)，令偏導(dǎo)數(shù)為 0,可得In Lnln f(Vj| )i 1L取極大值時(shí)，lnL也同時(shí)取極大值。求式(1.2 )(1.2)(1.3)解上式可得的極大似然估計(jì)ML。2系統(tǒng)參數(shù)的極大似然估計(jì)Newto n-Ra phs on 法實(shí)際上就是一種遞推算法，可以用于在線辨識。不過它是一種依 每L次觀測數(shù)據(jù)遞推一次的算法，現(xiàn)在我們討論的是每觀測一次數(shù)據(jù)就遞推計(jì)算一次參數(shù)估 計(jì)值得算法。本質(zhì)上說，

3、它只是一種近似的極大似然法。設(shè)系統(tǒng)的差分方程為式中a(z 1)y(k)b(z1)u(k)(k)(2.1)1a(z )b(z1)bo1aizb1Z 1n.anz.bnz n因?yàn)?k)是相關(guān)隨機(jī)向量，故(2.1)可寫成1a(z )y(k) b(z1)u(k)c(z 1) (k)(2.2 )式中(k)是均值為0a1, >a>b0>11 >C1 >設(shè)待估參數(shù)C(z1) (k)C(z 1)1 C1z 1的高斯分布白噪聲序列。Cn和序列ai并設(shè)y(k)的預(yù)測值為y(k)a1y(k 1)cie(k 1)(k)Cn多項(xiàng)式(k)的均方差anb0bnC1an y(kCne(k n)式

4、中e(k i)為預(yù)測誤差；ai， bi，Ci 為 ai，或者e(k)nCe(ki 1i)(1y(k) y(k)(1a1 z 1(C1 z 12C2 zy(k)1C1 zCn Z(2.3)(2.4) a(z 1)， b(z 1)和 C(z 1)中的系數(shù)都是未知參數(shù)。znTCn(2.5)n) b0u(k)bnu(k n)(2.6)bi，Ci的估值。naiy(ki 1anz n)y(k) (boCnz n)e(k)n)e(k) = (1 a1 z 11zi)預(yù)測誤差可表示為nbu(k i)i 0b1 zan Zbnzn)u(k)n)y(k)(2.7)bnz n)u(k)(2.8)1C(z )e(k)

5、a(:z1)y(k) b(z1)u(k)a(z1)1a-i z 1an zb(z1)bo bl z 1bn zC(z1)1C| z 1Cn z假定預(yù)測誤差e(k)服從均值為0的高斯分布，并設(shè)序列因此預(yù)測誤差e(k)式中nnne(k)具有相同的方差(bo b1滿足關(guān)系式(2.9)111為e(k)與c(z )， a(z )和b(z )有關(guān)，所以2是被估參數(shù) 的函數(shù)。為了書寫方便,把式(2.9)寫成c(z 1)e(k) a(z 1)y(k) b(z 1)u(k)n) b0u(k 1) biu(k 1)e(k) y(k) aiy(k 1)any(k或?qū)懗蒪nU(k n) cie(k1)Cn(kn),

6、k n 1,n2,(2.10 )(2.11 )e(k)ny(k)aiy(ki 1i)nbu(ki 0i)i令 k=n+1,n+2,cie(k i)1N個(gè)方程式寫成向量(2.12 )-矩陣形式式中eNYnNa1y(n1)e(n1)y(n2),eNe(n2)Jan boy(nN)e(nN)bnYn把這,n+N,可得e(k)的N個(gè)方程式,(2.13 )y(n)y(1)u(n1)u(1)e(n)e(1)y(n1)y(2)u(n2)u(2)e(n 1)e(2)y(nN 1)y(N)u(nN)u(N)e(n N 1)e(N)0的高斯噪聲序列，高斯噪聲序列的概率密度函數(shù)為因?yàn)橐鸭俣╡(k)是均值為1式中y為

7、觀測值,f (22和re農(nóng) Jy(y m)22)2 2(2.14 )m為y的方差和均值，那么1 1 2exp Le(k)(22)22對于e(k)符合高斯噪聲序列的極大似然函數(shù)為(2.15 )L(Yn , ) Le(n 1),e(n 2),e(n N)1 1 2 2 exp Le (n 1) e (n 2) (2 V 2fe(ne2(n1) fe(n 2)1N)N exp(2于fe(n n)22 eN eN )(2.16 )1 L(YN , ) exp(2 2)2對上式(2.17 )等號兩邊取對數(shù)得1N2尸lnL(YN ,)ln 一(2In e初或?qū)憺镮n L(YnNln22求 In L(Yn

8、,)對22的偏導(dǎo)數(shù), ln L(Yn| ,)令其等于2 1 nNk(k)式中2越小越好，因?yàn)楫?dāng)方差小)T(Yn0,2丄 N 2kn12 k n2最小時(shí)，(2.17 )T 、2 eNeN)Nln2(2.18 )1 n廠k可得Ne2(k)1Ne2(k)1 1e2(k)1(2.19 )(2.20 )(2.21 )n Ne2(k)1e2(k)最小，即殘差最小。因此希望(2.22 )2的估值取最2 2min JN因?yàn)槭?2.10 )可理解為預(yù)測模型，而 e(k)可看做預(yù)測誤差。因此使式(是使誤差的平方之和最小，即使對概率密度不作任何假設(shè)，這樣的準(zhǔn)則也是有意義的?？砂碕最小來求a1, ,a,b0, bn,

9、C1, cn的估計(jì)值。(2.23 )2.22 )最小就因此由于e(k)式參數(shù)a, ,a,b0, bn,c, cn的線性函數(shù)，因此J是這些參數(shù)的二次型函數(shù)。求使lnL(YN ,)最大的，等價(jià)于在式(2.10 )的約束條件下求使J為最小。由于J對Ci是非線性的，因而求 J的極小值問題并不好解，只能用迭代方法求解。求J極小值的常用迭代算法有拉格朗日乘子法和牛頓-拉卜森法。下面介紹牛頓-拉卜森法。整個(gè)迭代計(jì)算步驟如下:(1)確定初始的 0值。對于0中的a1, ,a,b0, bn可按模型e(k) a(z1 1)y(k) b(z )u(k)(2.24 )用最小二乘法來求，而對于0中的5' Cn可先

10、假定一些值。(2.29 )給出(2)計(jì)算預(yù)測誤差e(k) y(k) y(k)(2.25 )12kn Ne2(k)n 1并計(jì)算1 'NkN2e1(k)(2.26 )計(jì)算J的梯度和海賽矩陣式中e(k)e(k)ac1e(k1)同理可得e(k)1e(k)(2.27 )e(k)e(k2e(k)an boe(k) e(k)aibnCie(k)Cn尹(k) a1y(k 1)Cne(kn)y(ki)aiC2any(ke(k 2)ain)bou(k)biu(k 1)n)bnu(k n)(2.28 )e(k)aiy(ki)Cjj 1e(kaij)e(k)biu(ki)nCj +j 1b(2.30 )e(k

11、)Cie(ki)nCj1Cie(k j)(2.31 )將式（2.29 ）移項(xiàng)化簡，有y(k i)罟Cje(k j)ainCj0aie(k j)(2.32 )因?yàn)閑(kj)e(k)z j(2.33 )由e（k j）求偏導(dǎo)，故e(kaij)e(k)z jai(2.34 )將(2.34 )代入（2.32 ）,所以y(ki)nCjj 0aie(k j)Cjj 0e(k)z所以得同理可得（e(k)ai j(2.35 )c(z1)1CZCnZaiy(ki)(2.36 )2.30 )和(2.31 )為c")罟u(ki)(2.37 )e(ki)(2.38 )C(z1) ek (i j)ajyk (i

12、j)jy(k i)(2.39 )根據(jù)（2.36 ）構(gòu)造公式將其代入(2.36 ),可得1041丄(2.46 )0c(z1) ek (i j)ajc(z1)旦)ai(2.40 )消除c(z 1)可得e(k)e(k ij)e(k i1)aiaja11( 2.38 )式e(k)e(k ij)e(k i)bibjb0e(k)e(k ij)e(k i1)CiCjC1同理可得(2.37 )和(2.41 )(2.42 )(2.43 )式( 2.29 )、式( 2.30 )和式(2.31 )均為差分方程，這些差分方程的初始條件為可通過求解這些差分方程，分別求出e(k)關(guān)于a, ,a,b0, bn,G,G的全部

13、偏導(dǎo)數(shù)，而這些偏導(dǎo)數(shù)分別為 y(k)，u(k)和e(k)的線性函數(shù)。下面求關(guān)于的二階偏導(dǎo)數(shù)，即Tn N e(k) e(k)nN2e(k)科k n 1(2.44 )當(dāng)接近于真值時(shí),e(k)接近于0。在這種情況下，式(2.44 )等號右邊第2項(xiàng)接近于0,可近似表示為(2.45 )J n N e(k) e(k) T2k n 1則利用式(2.45 )計(jì)算比較簡單。(4)按牛頓-拉卜森計(jì)算的新估值1，有重復(fù)(2)至(4)的計(jì)算步驟，經(jīng)過r次迭代計(jì)算之后可得 r，近一步迭代計(jì)算可得如果10(2.47 )(2.48 )則可停止計(jì)算，否則繼續(xù)迭代計(jì)算。式(2.48 )表明，當(dāng)殘差方差的計(jì)算誤差小于O.01%時(shí)

14、就停止計(jì)算。這一方法即使在噪聲比較大的情況也能得到較好的估計(jì)值實(shí)驗(yàn)內(nèi)容設(shè)SISO系統(tǒng)差分方程為y(k) aiy(k 1)a2y(k 2) b1u(k 1) b2u(k 2)(k)其中極大似然估計(jì)法默認(rèn) e(k)為:e(k) C(z1) (k)(k) G (k1) C2 (k 2)辨識參數(shù)向量為=a1a2blb2c1c2t式中，(k)為噪聲方差各異的白噪聲或有色噪聲;u(k)和y(k)表示系統(tǒng)的輸入輸出變量。四實(shí)驗(yàn)流程圖Y停機(jī)1a.5 RS- *茯的程沖流耗ra五代碼實(shí)現(xiàn)程序如下:U=1.147 0.201 -0.787 -1.584 -1.052 0.866 1.152 1.573 0.626

15、0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 1.0991.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.8900.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.6090.090 -0.813 -0.428 -0.848 -0.410 0.048 -1.099 -1.108 0.259-1.627 -0.52

16、8 0.203 1.204 1.691 -1.235 -1.228 -1.267 0.3090.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.4630.786 1.161 0.850 -1.349

17、 -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.6810.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.2421.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

18、 0.206 1.196 1.0131.518 -0.553 -0.987 0.167 -1.445 0.630 1.255 0.311 -1.7260.975 1.718 1.360 1.667 1.111 1.018 0.078 -1.665 -0.7601.184 -0.614 0.994 -0.089 0.947 1.706 -0.395 1.222 -1.3510.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.20

19、40.926 -1.297 %輸入數(shù)據(jù)Y=0.086 2.210 0.486 -1.812 -3.705 -2.688 1.577 2.883 3.7051.642 0.805 -2.088 0.946 -0.039 1.984 -2.545 -1.727 -0.2312.440 3.583 2.915 1.443 3.598 0.702 2.638 3.611 -0.1681.732 0.666 2.377 -0.554 -2.088 2.698 0.189 -1.633 -2.0101.716 -1.641 -1.885 1.061 -0.968 2.911 3.088 -1.629 -1

20、.5333.030 0.614 -1.483 -1.029 -1.948 -1.066 -0.113 -2.144 -2.6260.134 -3.043 -1.341 0.338 2.702 3.813 -1.924 -2.813 -1.7953.002 1.027 1.027 2.755 3.584 1.737 -0.837 -0.617 1.7032.045 -2.886 -0.542 -2.991 -1.859 3.045 0.068 -0.375 0.4511.036 0.153 -0.474 2.512 -2.681 -0.954 -0.307 0.628 -0.270-0.277

21、0.983 -0.288 0.846 1.325 0.723 1.750 1.401 1.3400.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.5951.809 -2.204 -2.330 -0.454 1.290 2.080 -1.990 -0.770 1.240-0.252 3.137 -2.379 1.

22、206 1.221 -1.977 2.471 -1.680 1.1481.816 0.055 -1.856 0.269 -1.323 -2.486 1.958 0.823 2.4812.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.83

23、9 1.917 0.865 3.1803.261 -2.755 -0.536 -1.171 -0.905 -3.303 -0.834 2.490 3.0390.134 1.901% 輸岀數(shù)據(jù)n a=2; nb=1; nc=2;d=1;nn=max (na,n c);L=size(Y,2);yfk=zeros (nn ,1); %yf(k-i) ufk=zeros( nn ,1); %uf(k-i)參數(shù)估計(jì)初值xiefk=zeros (n c,1); %vf(k-i) thetae_1=zeros( na+n b+1+ nc,1); %P=eye( na+n b+1+ nc);for k=3:L

24、組建h (k)%構(gòu)造向量phi=-Y(k-1);-Y(k-2);U(k-1);U(k-2);xiek; %xie=Y(k)-phi'*thetae_1;phif=-yfk(1: na);ufk(d:d+nb);xiefk;%遞推極大似然參數(shù)估計(jì)算法K=P*p hif/(1+ phif* P*p hif);thetae(:,k)=thetae_1+K*xie;P=(eye( na+nb+1+ nc)-K* phif)* P;yf=Y(k)-thetae( na+n b+2: na+n b+1+ nc,k)'*yfk(1: nc); %yf(k)uf=U(k)-thetae( na

25、+nb+2:n a+nb+1+ nc,k)'*ufk(1: nc); %uf(k)xief=xie-thetae( na+n b+2:na+n b+1+ nc,k)'*xiefk(1: nc); %vf(k)%更新數(shù)據(jù)thetae_1=thetae(:,k);for i=n c:-1:2xiek(i)=xiek(i-1);xiefk(i)=xiefk(i-1);endxiek(1)=xie;xiefk(1)=xief;for i=nn :-1:2yfk(i)=yfk(i-1);ufk(i)=ufk(i-1);endyfk(1)=yf;ufk(1)=uf;end thetae_1 figure(1) plot(1:L,thetae(1: na,:);xlabel('k'); ylabel('參數(shù)估計(jì) a');lege nd('a_1','a_2'); axis(0 L -2 2);figure(2) plot(1:L,thetae (n a+1: na