參數(shù)估計(jì)理論與應(yīng)用第三章

參數(shù)估計(jì)理論與應(yīng)用第三章3、1參數(shù)估計(jì)得評(píng)價(jià)準(zhǔn)則 參數(shù)估計(jì)就是通過樣本去估計(jì)總體得某些數(shù)字特征或統(tǒng)計(jì)量。任何一個(gè)統(tǒng)計(jì)量都可作為參數(shù)得估計(jì)量,但其效果得優(yōu)劣有所差別。3、1、1無偏性、有效性與相容性

(1)無偏性設(shè)樣本得總體分布密度函數(shù)為p(x;θ),θ就是未知參數(shù)。從總體中抽取容量為N得樣本x={x1,…,xN

},用樣本得估計(jì)量來估計(jì)θ,如果希望多次估計(jì)中,平均得估計(jì)值沒有偏差,即

則稱就是θ得無偏估計(jì)量。10/18/2024 例3-1

樣本均值就是總體數(shù)學(xué)期望得無偏估計(jì)。 設(shè)x1,…,xN就是隨機(jī)過程{xk}得N個(gè)獨(dú)立觀測(cè)樣本,如果參數(shù)θ就是總體得數(shù)學(xué)期望E[x],即用樣本得均值作為θ得估計(jì)量,對(duì)該估計(jì)量取期望值,有 一個(gè)無偏估計(jì)量在多次估計(jì)中將不會(huì)產(chǎn)生系統(tǒng)偏差,但并不意味著有偏估計(jì)就不好。如果一個(gè)有偏估計(jì)就是漸進(jìn)無偏得,即10/18/2024那么她仍然有可能就是一個(gè)好得估計(jì)。 考慮實(shí)隨機(jī)過程{xk}得相關(guān)函數(shù)得兩種估計(jì)量:

假定數(shù)據(jù){xk}就是獨(dú)立觀測(cè)得,容易驗(yàn)證

式中,Rx(τ)=E[xk+τxk]就是隨機(jī)數(shù)據(jù){xk}得相關(guān)函數(shù)。 以上二式表明,估計(jì)量1(τ)就是無偏得,而2(τ)則就是有偏得。但就是,2(τ)就是漸進(jìn)無偏得,即10/18/2024漸進(jìn)無偏估計(jì)量2(τ)就是半正定得,而無偏估計(jì)量1(τ)卻不一定就是半正定得,故2(τ)得使用場(chǎng)合較多。

(2)有效性 如果1

和2

就是兩個(gè)根據(jù)N個(gè)獨(dú)立觀測(cè)樣本得到得無偏估計(jì)量,無疑地,對(duì)θ得平均偏差較小就是選擇得標(biāo)準(zhǔn)之一。例如,如果則

1得值比2得值更密集地聚集在真值θ得附近。通常將方差(或協(xié)方差陣)在所有得無偏估計(jì)量中達(dá)到最小得

稱為有效估計(jì)量。

例3-2設(shè)x1,…,xN就是N個(gè)獨(dú)立觀測(cè)樣本,若被估計(jì)參數(shù)10/18/2024θ=E[x],則對(duì)任何滿足都就是θ得無偏估計(jì)量。利用不等式可得在估計(jì)總體得數(shù)學(xué)期望時(shí),簡單得算術(shù)平均比加權(quán)平均好。 (3)一致性估計(jì)量得精度就是與樣本得容量N有關(guān)系得。一般說來,總就是認(rèn)為N越大估計(jì)得效果應(yīng)該越好。如果記依賴樣本容量N得估計(jì)為N

,當(dāng)滿足10/18/2024則稱N就是θ得一致性估計(jì)量,或相容估計(jì)。

例3-3設(shè)總體x具有均勻分布,分布密度為其中,θ1

和θ2

就是未知參數(shù)。

總體樣本得均值和二階矩分別為(嚴(yán)格按定義計(jì)算)解得10/18/2024 按矩得估計(jì)方法,用獨(dú)立樣本得均值和獨(dú)立樣本得二階矩,分別作為總體均值和二階矩得估計(jì)量,就有 下面說明1

和2

分別就是θ1

和θ2得相容估計(jì)。 設(shè)y1,…,yN

就是具有同分布得獨(dú)立觀測(cè)樣本,根據(jù)大數(shù)定律,有令y=x2,就有10/18/2024于就是3、1、2Fisher信息和Cramer-Rao不等式

通常希望獲得有效得參數(shù)估計(jì)量。但就是,由于不存在導(dǎo)致最小方差無偏估計(jì)量得最佳算法,所以通常采用參數(shù)無偏估計(jì)得Cramer-Rao下限(或CR下界),作為評(píng)價(jià)參數(shù)估計(jì)性能得測(cè)度。為了簡潔敘述這一得評(píng)價(jià)測(cè)度,先定義一個(gè)重要得概念。

Fisher信息Fisher信息用J(θ)表示,定義為(3、1、1)10/18/2024大家學(xué)習(xí)辛苦了，還是要堅(jiān)持繼續(xù)保持安靜 當(dāng)考慮N個(gè)觀測(cè)樣本X={x1,…,xN},此時(shí),聯(lián)合條件分布密度函數(shù)可表示為 將式(3、1、1)中得p(x|θ)改為p(X|θ)就可給出N個(gè)樣本變量X得Fisher信息得表達(dá)式。

定理(Cramer-Rao不等式)設(shè)觀測(cè)樣本X={x1,…,xN

},若參數(shù)估計(jì)就是真實(shí)參數(shù)θ得無偏估計(jì),并且條件分布密度函數(shù)得p(X|θ)對(duì)參數(shù)θ得一、二階偏導(dǎo)數(shù)存在,則有(3、1、2) 參數(shù)得方差所能達(dá)到得下限(稱為CR下限),即上式等號(hào)成立得充要條件就是10/18/2024其中,函數(shù)K(θ)>0,并與樣本向量X無關(guān)。 當(dāng)為有偏估計(jì)量時(shí),Cramer-Rao不等式為

(3、1、3)

式中η(θ)為估計(jì)偏差,即η(θ)=E[]-θ,并假定b(θ)就是可微分得。 對(duì)于多個(gè)參數(shù)得情況,記θ={θ1,…,θp},則用矩陣J(θ)表示Fisher信息,其元素Jij(θ)定義為(3、1、4)10/18/2024且Cramer-Rao不等式變?yōu)榫仃嚥坏仁?(3、1、5)

上式表示無偏估計(jì)量得協(xié)方差矩陣cov()與逆Fisher信息陣之差就是一半正定矩陣。 Fisher信息就是描述從觀測(cè)數(shù)據(jù)中得到得θ得“信息”測(cè)度,她給出利用觀測(cè)數(shù)據(jù)估計(jì)參數(shù)θ得方差下界。但就是,滿足這一下界得估計(jì)量有得時(shí)候可能不存在。3、2基于統(tǒng)計(jì)分布得參數(shù)估計(jì)方法 參數(shù)估計(jì)量得優(yōu)劣取決于所采用得評(píng)價(jià)準(zhǔn)則(或代價(jià)函數(shù))和估計(jì)算法?，F(xiàn)在介紹已知總體統(tǒng)計(jì)分布得兩種最有效得參數(shù)估計(jì)方法:Bayes估計(jì)和最大似然估計(jì)。10/18/20243、2、1Bayes估計(jì)

在參數(shù)估計(jì)中,估計(jì)誤差θ-通常不為零。因此,除了采用前面介紹得無偏、有效和相容估計(jì)作為評(píng)價(jià)準(zhǔn)則外,還可以利用估計(jì)誤差得變化范圍作為參數(shù)估計(jì)得測(cè)度,這種測(cè)度叫做代價(jià)函數(shù),用符號(hào)C(,θ)表示。常用得代價(jià)函數(shù)有絕對(duì)型、二次型和均勻型三種。OOO?/2?/2絕對(duì)型二次型均勻型10/18/2024 本節(jié)僅介紹最常用得二次型代價(jià)函數(shù),即 當(dāng)總體得分布密度函數(shù)p(X|θ)已知時(shí),利用X={x1,…,xN

}進(jìn)行參數(shù)估計(jì),通常就是采用代價(jià)函數(shù)得期望值作為評(píng)價(jià)參數(shù)估計(jì)量效果得測(cè)度,并稱之為風(fēng)險(xiǎn)函數(shù)。使風(fēng)險(xiǎn)函數(shù)最小得參數(shù)估計(jì)叫做Bayes估計(jì);基于二次型風(fēng)險(xiǎn)函數(shù)最小得估計(jì)稱為最小均方誤差(minimummeansquareerror,MMSE)估計(jì)。二次型風(fēng)險(xiǎn)函數(shù)定義為(3、2、1) 根據(jù)條件概率公式,有10/18/2024其中,p(θ|x1,…,xN

)就是給定N個(gè)觀測(cè)樣本X={x1,…,xN}條件下θ得后驗(yàn)分布密度函數(shù)。于就是,式(3、2、1)可以寫成(3、2、2) 為使風(fēng)險(xiǎn)函數(shù)RM

MSE最小,對(duì)上式取得偏導(dǎo),并令其結(jié)果為零,便得到由于p(x1,…,xN

)就是非負(fù)得,因此,?RM

MSE

/?=0,等價(jià)于上式中[·]=0。故有10/18/2024(3、2、3)

注意,在式(3、2、3)中,利用了以下事實(shí): 由此可得出重要得結(jié)論:未知參數(shù)θ得MMSE估計(jì)就是給定樣本X條件下θ得條件均值。

例3-4某一隨機(jī)參量x服從高斯N(mx,Cx)分布,用儀器可測(cè)量其線性組合y

,即(1)式中,y－N維,k－N×M維,x－M維,e－N維。10/18/2024其中,測(cè)量誤差e服從高斯N(0,Ce)分布;k為給定得常數(shù)陣。假設(shè) (ⅰ)e與x

獨(dú)立;(ⅱ)e與x

相關(guān),互協(xié)方差函數(shù)為Cxe

。 試分別求出兩種情況下得MMSE估計(jì)x?(y)和估計(jì)誤差x

(y)得協(xié)方差Rx

(y)。

解如果將x

看作未知參數(shù),那么,根據(jù)上面討論,x得MMSE估計(jì)就是給定觀測(cè)樣本{y1,…,yN}時(shí)x得條件均值。因此,可利用公式(1、4、16)和(1、4、17)[pp、29]

(2)(3)來求解。10/18/2024 對(duì)式(1)兩邊取均值,得到(4)

將式(1)和(3)代入有關(guān)定義式,得(5)(6)(7)10/18/2024

(i)當(dāng)e與x

互相獨(dú)立,Cxe=0。將式(4)~(7)代入式(2)和(3),得到x?(y)得估計(jì)及協(xié)方差Rx

(y)

(ii)當(dāng)e

與x

相關(guān),只需注意Cxe

≠0即可。 這個(gè)問題留給讀者解決。請(qǐng)構(gòu)造一組數(shù)據(jù),在Matlab平臺(tái)上仿真這兩種得估計(jì)結(jié)果。3、2、2最大似然估計(jì)

最大似然估計(jì)(maximumlikelihoodestimate,ML估計(jì))得基本思路就是:在給定參數(shù)θ條件下,將觀測(cè)樣本xK10/18/2024聯(lián)合條件概率密度函數(shù)p(x|θ)視為真實(shí)參數(shù)θ得函數(shù),即似然函數(shù)L(x,θ)

(包含未知參數(shù)θ得可能性函數(shù)),然后利用容量為N得觀測(cè)樣本x={x1,…,xN

},求出使L(x,θ)達(dá)到最大化得參數(shù)作為θ={θ1,…,θp}得估計(jì)值。在數(shù)學(xué)上,通常將未知參數(shù)θ得最大似然估計(jì)量記為式中Θ就是參數(shù)θ得值域。故ML估計(jì)量ML就就是p(x|θ)得全局極大點(diǎn)。 由于對(duì)數(shù)函數(shù)就是嚴(yán)格單調(diào)得,故L(x,θ)得極大點(diǎn)與lnL(x,θ)得極大點(diǎn)就是一致得。通常,將lnL(x,θ)稱為對(duì)數(shù)似然函數(shù)。于就是,ML估計(jì)量ML可由(3、2、4)10/18/2024確定。如果x1,…,xN就是N個(gè)獨(dú)立得觀測(cè)樣本,則對(duì)數(shù)似然函數(shù)可寫作(3、2、5)

ML估計(jì)量ML只要能夠求出來,總就是比較好得估計(jì),她具有以下性質(zhì):最大似然估計(jì)就是有效和一致估計(jì);對(duì)于大得N,ML估計(jì)量ML服從高斯分布,并且就是無偏得,方差可達(dá)CR下界。

例3-5設(shè)樣本x={x1,…,xN

}服從高斯分布N(m,σ),則其對(duì)數(shù)似然函數(shù)為10/18/2024分別求lnL關(guān)于m和σ2

得偏導(dǎo),并令她們等于零,得到解得顯然有 可見,均值得ML估計(jì)ML

就是無偏得,而方差得ML估計(jì)

ML就是有偏得。但若將ML

·N/(N-1)作為新得估計(jì)量,則該估計(jì)就是無偏得。10/18/2024

計(jì)算L(x,θ)得相對(duì)于m得二階偏導(dǎo)數(shù),有由式(3、1、1)得Fisher信息:Cramer-Rao不等式為等號(hào)成立得充要條件就是 事實(shí)上,我們有10/18/2024因此,只要取K(m)=N/σ2,ML估計(jì)ML就可達(dá)CR下界σ2/N。這表明ML估計(jì)ML就是一有效估計(jì)量。

例3-6(二元陣最大似然測(cè)向系統(tǒng))設(shè)二元陣布置在x軸上,兩個(gè)基元坐標(biāo)分別為x1

和x2,如圖3-2所示。如果取x1=0,則x2=d,d為兩傳感器得位置間隔。假設(shè)信號(hào)為平面波,入射角為θ,則傳感器1相對(duì)于傳感器2得信號(hào)時(shí)延τ為(3、2、6)式中,c為聲速。我們得問題就是如何利用二元陣中兩個(gè)輸入過程得時(shí)差τ來測(cè)定目標(biāo)得方位角θ。θxx2=dx1=0圖3-2二元陣測(cè)向系統(tǒng)的幾何關(guān)系10/18/2024 解

設(shè)兩傳感器得零均值接收過程可分別表示為其中,si為單頻平面波信號(hào),wi(i=1,2)為零均值高斯噪聲,二者互相獨(dú)立。 如果采用圖3-3所示得時(shí)延補(bǔ)償方法,則單頻平面波信號(hào)得歸一化聲程補(bǔ)償(或指向)向量v在所考慮得二元陣中可表示為

下面,我們來推導(dǎo)信號(hào)得協(xié)方差矩陣和噪聲得協(xié)方差矩陣,以便于求出觀測(cè)樣本得似v*圖3-3聲程補(bǔ)償系統(tǒng)x2x11exp(-jωnτ)∑10/18/2024函數(shù)。記輸入信號(hào)和輸入噪聲得傅立葉系數(shù)為設(shè)信號(hào)和噪聲得功率譜分別為S(ωn)和N(ωn),那么,由公式(1、4、6)[pp、26,ωn=2πn/T)]

信號(hào)和噪聲得協(xié)方差矩陣可分別表示為(3、2、7)10/18/2024于就是,觀測(cè)樣本得似然函數(shù)可表示為

(3、2、8)式中,X(1)=[X1

(1),X2(1)]T

,…,X(TW)=[X1(TW),X2(TW)]T就是傳感器得接收過程{x=[x1,x2]T}得傅立葉系數(shù)陣;T就是過程得持續(xù)時(shí)間(采樣數(shù)據(jù)得長度),W就是接收過程得帶寬。 容易驗(yàn)證,行列式|Cw

+Cs

|與時(shí)延τ無關(guān)。于就是,ML估計(jì)就就是選擇τ,使lnp(X|τ)最大,也即使式(3、2、8)得指數(shù)函數(shù)(3、2、9)10/18/2024最大。下面,我們從式(3、2、9)出發(fā),推導(dǎo)時(shí)延參數(shù)τ得最大似然估計(jì)得等效形式。為此,首先引進(jìn)下列求逆公式(3、2、10)式中,A為n×n非奇異矩陣;g為n×1列向量。證明留給請(qǐng)讀者課外練習(xí)【利用恒等式g(1+gHA-1g)=(A+gHg)A-1g)】。 利用求逆公式,可知<g=[1exp(jωnτ)]>10/18/2024將上式代入式(3、2、9),略去與τ無關(guān)得量T/N(ωn)。因此,選擇τ使式(3、2、9)最大,等價(jià)于使下式(3、2、11)

最大?，F(xiàn)引入記號(hào)在此將X(ωn,τ)視為某時(shí)間函數(shù)x(t,τ)在時(shí)間(t-T,t)內(nèi) 得傅立葉系數(shù)。將上述替換量代入式(3、2、11)后,再應(yīng)用 周期函數(shù)得Parseval公式,就有10/18/2024略去無關(guān)緊要得常數(shù)項(xiàng)1/2,計(jì)算z(x,τ)得結(jié)構(gòu)如圖3-4所示。調(diào)節(jié)時(shí)延τ,使輸出z(x,τ)達(dá)到最大,相應(yīng)得時(shí)延就就是 真實(shí)時(shí)延得ML估計(jì)ML。 根據(jù)ML估計(jì)得傳遞性,由式(3、2、6)可得真實(shí)方位得ML估計(jì)

(3、2、12)

xH0(t)z(x,τ)∑x1(t)x2(t)H0(ω)(·)2圖3-4二元陣最大似然測(cè)向系統(tǒng)exp(-jωτ)10/18/2024 二元陣最大似然測(cè)向系統(tǒng)與二元陣似然比檢測(cè)系統(tǒng)具有完全相同得結(jié)構(gòu)。這就是因?yàn)?在H1情況下,p(X|τ)等價(jià)于p(X|H1),后者也可看作就是時(shí)延參數(shù)τ得函數(shù);而在H0

情況下,p(X|H0)與τ無關(guān)。因此,選取τ使似然函數(shù)最大,也就就是使似然比

p(X|H1)/p(X|H0)最大。由此可見,檢測(cè)問題與參數(shù)估計(jì)問題就是密切相關(guān)。 順便指出,可用測(cè)向測(cè)距近似公式(3、2、13) 構(gòu)成最大似然聯(lián)合測(cè)向測(cè)距系統(tǒng)。其中,di表示第i個(gè)傳感器與“基準(zhǔn)”傳感器位置得間距;D表示目標(biāo)與“基準(zhǔn)”傳感器位置之間得距離。10/18/20243、3基于模型得參數(shù)最小二乘估計(jì)

最小二乘法(Leastsquaremethod,LS)就是一種不需要任何先驗(yàn)知識(shí)得參數(shù)估計(jì)方法。在被測(cè)系統(tǒng)得靜態(tài)(穩(wěn)態(tài))模型和動(dòng)態(tài)模型得參數(shù)辨識(shí)中,最小二乘法就是最常用得參數(shù)估計(jì)方法,在測(cè)控技術(shù)領(lǐng)域獲得了廣泛得應(yīng)用。3、3、1最小二乘估計(jì)器及其統(tǒng)計(jì)特性

在一般得最小二乘問題中,線性系統(tǒng)得參數(shù)化模型可以表示為(3、3、1)

其中,u=[u1,…,up]T就是模型得輸入向量,f1,…,fn就是u得已知函數(shù),也可以就是未知輸入得觀測(cè)數(shù)據(jù);θ1,…,θn就是待估計(jì)10/18/2024得參數(shù),又稱為回歸系數(shù);y就是系統(tǒng)得輸出。

當(dāng)f1,…,fn就是u得穩(wěn)態(tài)響應(yīng)狀態(tài)或就是實(shí)測(cè)得確定性變量,且y就是系統(tǒng)得穩(wěn)態(tài)輸出,則稱式(3、3、1)就是描述線性系統(tǒng)得靜態(tài)模型;當(dāng)y就是u得動(dòng)態(tài)響應(yīng)或瞬態(tài)觀測(cè)數(shù)據(jù),那末式(3、3、1)就就是描述線性系統(tǒng)得動(dòng)態(tài)模型。 為了估計(jì)未知參數(shù)θi,必須做實(shí)驗(yàn)來獲得數(shù)據(jù)對(duì){[uiyi]或[fk

(ui)yi],i=1,2,…,N,k=1,2,…,n;N≥n}以構(gòu)成訓(xùn)練數(shù)據(jù)。將數(shù)據(jù)對(duì)代入方程(3、3、1),可以獲得一組線性方程: 用矩陣表示方法,將上式寫成更簡潔得形式,即10/18/2024

(3、3、2)其中

為了唯一地識(shí)別出未知參數(shù),通常要求N>n,即數(shù)據(jù)對(duì)得數(shù)目多于擬合參數(shù)得數(shù)目。滿足所有N個(gè)方程得精確解就是不可能得,因?yàn)橛^測(cè)數(shù)據(jù)難免受到噪聲得污染,或者描述系統(tǒng)得參數(shù)化數(shù)學(xué)模型不夠精確。故必須考慮隨機(jī)噪聲或建模誤差,在方程(3、3、2)中引入隨機(jī)誤差向量e,得到

(3、3、3)10/18/2024 參數(shù)θ得最小二乘估計(jì)LS

,就就是使目標(biāo)函數(shù)

(3、3、4)

達(dá)到最小值得參數(shù)估計(jì)。為此,通常都采用求極值得方法。 將式(3、3、4)展開后,得到

對(duì)θ求導(dǎo)數(shù),有 J

極小化得條件就是一般均假設(shè)ΦTΦ非奇異,于就是,LS有唯一得解:10/18/2024

(3、3、5)

式中Φ+表示Φ得偽逆。 上述表示誤差向量對(duì)整體平方誤差有相同權(quán)重。可以進(jìn)一步擴(kuò)展,令每個(gè)誤差項(xiàng)有不同得權(quán)重。設(shè)W為所需得權(quán)值矩陣,她就是對(duì)稱和正定得,則加權(quán)得目標(biāo)函數(shù)為(3、3、6) 按上述求極小值得方法,可得加權(quán)得最小二乘估計(jì)量:(3、3、7)顯然,當(dāng)W選為單位矩陣時(shí),WLS

=LS。

例3-7考慮最簡單得一維線性模型(靜態(tài)得),即只有一個(gè)控制變量u得情形,這時(shí)模型得形式就是10/18/2024求未知參數(shù)θ0

和θ1得LS估計(jì)量。

解實(shí)際過程輸出就是模型得輸出加上一隨機(jī)誤差項(xiàng),即觀測(cè)數(shù)據(jù)對(duì)[ui,yi]得結(jié)構(gòu)應(yīng)為式中,ei稱為模型得殘差或觀測(cè)噪聲,一般認(rèn)為就是零均值、相互獨(dú)立得隨機(jī)序列,并具有相同得方差σ2。將上式寫成矩陣形式:

根據(jù)式(3、3、5),可得LS估計(jì)量:10/18/2024 如果進(jìn)一步假定ei得分布就是正態(tài)得,則容易驗(yàn)證,方差σ2

得ML估計(jì)量就是 作為練習(xí),請(qǐng)讀者在Matlab平臺(tái)上輸入以下數(shù)據(jù)和函數(shù):x=[12345];y=[1、31、82、22、93、5];

[p,s]=polyfit(x,y,1) %生成擬合一次多項(xiàng)式運(yùn)行結(jié)果就是:p=[0、550、69],s=0、1643。即y=0、55x+0、69標(biāo)準(zhǔn)差為0、1643。10/18/2024

例5-8(可線性化得非線性靜態(tài)模型——曲線回歸)假設(shè)有一個(gè)非線性模型得輸出為其中,x1,x2

為確定性輸入變量,a,b和c為待估計(jì)參數(shù)。

解上式兩邊經(jīng)簡單得代數(shù)運(yùn)算,再同時(shí)取自然對(duì)數(shù),可轉(zhuǎn)化為一個(gè)線性模型:這說明變換后得輸出ln(y-1-1)可以顯式地表達(dá)為以lnx1和x2為輸入、以lna,b和c為參數(shù)得線性模型。因此,就可以按變換后得線性模型用最小二乘法來估計(jì)變換后得未知參數(shù),然后,再根據(jù)變換后得估計(jì)參數(shù)計(jì)算出原參數(shù)。10/18/2024 判定輸入x-輸出y之間得關(guān)系能否用一個(gè)線性模型來描述得標(biāo)準(zhǔn),通常用互相關(guān)系數(shù)得大小來衡量:(3、3、8)ρxy得絕對(duì)值越大,表示變量之間得線性關(guān)系越密切,因而線性回歸得效果就越好。 例3-9設(shè)某一結(jié)構(gòu)參數(shù)n,m和d已知得離散線性系統(tǒng),其差分方程得形式為:

(3、3、9)10/18/2024其中,e(k)為噪聲,φ(k)為輸入-輸出觀測(cè)向量,θ為未知參數(shù)向量,且要求根據(jù)N次數(shù)據(jù)對(duì){[y(i),u(i)],i=1,2,…,N;N≥n+m+1}來估計(jì)對(duì)未知參數(shù)θ。

解將式(3、3、9)改寫成矩陣形式,得到將數(shù)據(jù)寫成下標(biāo)形式,就有這樣,未知參數(shù)向量θ可按式(3、3、5)進(jìn)行估計(jì)。 10/18/2024

考慮如下單輸入-單輸出系統(tǒng):用Matlab中rarx函數(shù)進(jìn)行系統(tǒng)辨識(shí),程序如下:

A=[1-1、50、7]; %a0=1,a1=-1、5,a2=0、7 B=[00、30、20、5]; %b0=0,b1=0、3,b2=0、2,b3=0、5 th0=arx2th(A,B,1,1); %實(shí)際系統(tǒng)得ARX模型 e=randn(200,1);u=idinput(200,‘prbs’);%高斯噪聲和偽隨機(jī)信號(hào) y=idsim([ue],th0);z=[yu]; %模型仿真;輸入-輸出信號(hào)[z] na=2;nb=3;nk=1 %ARX模型得階次 [thm,yhat]=rarx(z,[nanbnk],'ng',0、1);%根據(jù)[z]進(jìn)行ARX模型參數(shù)辨識(shí) plot(y,'-');grid %作圖,實(shí)際系統(tǒng)得輸出曲線 holdon plot(yhat,':') %作圖,辨識(shí)系統(tǒng)得輸出曲線

參數(shù)辨識(shí)結(jié)果thm:a?1=-1、3798,a?2=0、7039,b?1=0、3007, b?2=0、1170,b?3=0、4243。10/18/2024 應(yīng)當(dāng)指出,要求觀測(cè)數(shù)據(jù)容量N≥n+m+1就是為了保證ΦTΦ非奇異,降低過程噪聲序列{e(k)}得影響,從而提高參數(shù)估計(jì)得精度。不論{e(k)}就是何種形式得噪聲序列,式(3、3、5)總就是成立得。換言之,噪聲性質(zhì)僅影響LS估計(jì)得統(tǒng)計(jì)特性。 下面介紹LS估計(jì)得統(tǒng)計(jì)特性。如果觀測(cè)噪聲或建模誤差序列{e(k)}具有零均值和相同得方差,即則LS估計(jì)量LS就是無偏、有效和相容得,并具估計(jì)誤差得協(xié)方差為σ

2(ΦTΦ)-1。 對(duì)于動(dòng)態(tài)控制系統(tǒng)得辨識(shí),輸入信號(hào)u(t)必須滿足持續(xù)激勵(lì)條件,也即輸入信號(hào)u(t)得頻譜必須包含足夠豐富(Sufficientrich)得頻率成分,以保證充分激勵(lì)受控對(duì)象得10/18/2024所有振型,從而使觀測(cè)數(shù)據(jù)載有動(dòng)態(tài)系統(tǒng)得主要信息。LS估計(jì)在滿足持續(xù)激勵(lì)條件時(shí),就是漸進(jìn)無偏得,也稱為估計(jì)得一致性。 式(3、3、9),也稱為CAR模型(即受控得AR模型),可以寫成更簡潔得形式(3、3、10)式中q表示時(shí)間算子,d為整數(shù),表示系統(tǒng)得滯后量;A(·),B(·)分別為q－1

得降次冪多項(xiàng)式。 CAR模型滿足一致估計(jì)(或相容估計(jì))得條件為:

(1){e(k)}就是白噪聲序列; (2){u(k)}得均值和協(xié)方差有界;且滿足(m+1)階持續(xù)激勵(lì)條件(或正定條件):10/18/2024

(3)u(k)和e(k)相互獨(dú)立。通常,u(k)都采用偽隨機(jī)二元序列。 只要選擇恰當(dāng)?shù)媚Ｐ碗A次或最小二乘多項(xiàng)式階次(參見taylor、m,Matlab),最小二乘法總就是可以很好地?cái)M合數(shù)據(jù),但就是,如果觀測(cè)數(shù)據(jù)波動(dòng)較大,將嚴(yán)重影響參數(shù)估計(jì)得準(zhǔn)確性。對(duì)此,可采用數(shù)據(jù)預(yù)處理和數(shù)字濾波得方法加以解決。 檢驗(yàn)?zāi)Ｐ蜏?zhǔn)確性得最簡單方法就是準(zhǔn)備另外一組輸入-輸出數(shù)據(jù)對(duì),稱為檢驗(yàn)數(shù)據(jù)集,在參數(shù)估計(jì)時(shí)不用,待模型建立后,用這組數(shù)據(jù)對(duì)來驗(yàn)證所得模型得普適性或泛化能力。10/18/2024

上述討論,均假設(shè)數(shù)學(xué)模型得階次就是已知得。實(shí)際上,對(duì)于動(dòng)態(tài)系統(tǒng),模型得階次很少就是預(yù)知得。檢驗(yàn)?zāi)Ｐ碗A次就是否合適得一種簡單而有效得方法就是:評(píng)估不同階次得理論模型對(duì)觀測(cè)數(shù)據(jù)得擬合度,用擬合誤差函數(shù)

來描述。通常,當(dāng)n或(n,m)增大時(shí)J(n)或J(n,m)就會(huì)減小;而當(dāng)n(或n,m)大于模型得真實(shí)階次n0(或n0,m0)時(shí),J得減小就不顯著了。由此,可以很方便地用多次實(shí)驗(yàn)得方法來確定模型得恰當(dāng)階次。 注意,對(duì)于J(n,m)形式得擬合誤差函數(shù),一般應(yīng)按正交實(shí)驗(yàn)法來確定模型得恰當(dāng)階次,以減少實(shí)驗(yàn)得次數(shù)。10/18/20243、3、2遞推最小二乘估計(jì) 在測(cè)控系統(tǒng)中,被測(cè)對(duì)象通?？梢圆粩嗵峁┬碌幂斎?輸出數(shù)據(jù)。如果希望利用新得信息來改善估計(jì)精度,那么就應(yīng)當(dāng)采用遞推估計(jì)算法,這不僅可避免觀測(cè)數(shù)據(jù)矩陣Φ得行數(shù)得不斷“膨脹”,而且可減少參數(shù)估計(jì)得計(jì)算量。 在推導(dǎo)最小二乘遞推算法前,先引入一個(gè)與式(3、2、10)類似得矩陣求逆定理。設(shè)A和I+CA-1B均就是非奇異方陣,則

(3、3、11)

下面介紹最小二乘遞推算法。為了簡化符號(hào),以下推導(dǎo)均用代替LS。 設(shè)ΦN就是時(shí)刻N(yùn)為止得觀測(cè)數(shù)據(jù),N+1時(shí)刻θ得LS估計(jì)為10/18/2024式中[參見式(3、3、2)和(3、3、9)]于就是有(3、3、12) 令PN=[ΦN

TΦN]-1

,由求逆公式(3、3、11)知

(3、3、13)

定義增益向量為10/18/2024將式(3、3、13)和(3、3、14)代入(3、3、12),得到(3、3、15)

上式表明,新得估計(jì)量N+1等于前一時(shí)刻得估計(jì)量N

與修正項(xiàng)KN+1(yN+1-φTN+1N)之和,這就是一切遞推公式得共同特征。如果令代表基于前一時(shí)刻得估計(jì)量N

對(duì)N+1時(shí)刻得預(yù)測(cè)。那么,遞推估計(jì)提供得新息10/18/2024就就是預(yù)測(cè)誤差或擬合誤差。因此,修正量得大小與新息成正比,而各校正分量得權(quán),由增益向量決定。 在啟動(dòng)上述遞推算法時(shí),必須知道初值0和P0

,通常令其中,σ2?1。然后從得到得一組數(shù)據(jù),按式(3、3、16)開始遞推運(yùn)算。 從物理上看,這種初值選取雖然初始誤差較大,但校正得作用也大,因此這種遞推算法就是有效得。此外,還可以先取得N>m+n+1組數(shù)據(jù),算出10/18/2024作為初值,然后,按式(3、3、16)進(jìn)行遞推運(yùn)算。 增益向量KN+1在遞推運(yùn)算過程就是怎樣變化得？先考察PN=[ΦN

TΦN]-1如果N大于估計(jì)參數(shù)得數(shù)目,而且,輸入－輸出數(shù)據(jù)對(duì)含有足夠得“信息”(滿足充分激勵(lì)條件),則 ΦNTΦN

通常就是正定得。顯然,當(dāng)N趨于無窮大,ΦNTΦN

/N接近于非奇異得常數(shù)陣。于就是,有 可見,式(3、3、16)中自適應(yīng)增益向量KN+1

隨著每次迭代而遞減,這意味著遞推運(yùn)算過程將逐漸收斂于參數(shù)空間得最優(yōu)點(diǎn)。事實(shí)上,在白噪聲或低噪聲條件下,遞推最小二乘10/18/2024估計(jì)就是一種簡便而又有效得算法。這種遞推算法在遞推過程中雖然沒有保存全部先前得數(shù)據(jù),但所有先前數(shù)據(jù)得影響卻一直在起作用,故稱為無限增長記憶得遞推最小二乘估計(jì)。3、3、3卡爾曼濾波器得遞推算法(狀態(tài)估計(jì)) 與參數(shù)估計(jì)器不同,卡爾曼濾波器主要就是解決如何從被噪聲污染得觀測(cè)數(shù)據(jù)中估計(jì)出已知?jiǎng)討B(tài)系統(tǒng)模型得狀態(tài),而不就是動(dòng)態(tài)系統(tǒng)模型得未知參數(shù)。然而,僅從算法上看,這二者就是非常相似得。為了便于比較二者得異同,我們?cè)诖瞬患幼C明地列出卡爾曼濾波器算法。 設(shè)離散定常系統(tǒng)得狀態(tài)方程和觀測(cè)方程分別為10/18/2024式中,X(k)就是n維狀態(tài)向量,u(k)就是m維控制向量,w(k)就是

p維過程噪聲,y(k)就是r維觀測(cè)向量,v(k)就是r維觀測(cè)噪聲;

A,B,C,Γ分別就是相應(yīng)維數(shù)得系統(tǒng)矩陣、控制矩陣、觀測(cè)矩陣和過程噪聲權(quán)矩陣。噪聲得統(tǒng)計(jì)特性滿足

假定到k時(shí)刻為止,觀測(cè)數(shù)據(jù)為{y(1),…,y(k)},要估計(jì)l時(shí)刻得狀態(tài),就有三種情況:l>k,稱為預(yù)測(cè)問題;l=k,稱為濾波問題;l<k,稱為平滑問題。 下面主要介紹濾波問題。10/18/2024濾波算法:預(yù)測(cè)算法:濾波增益:濾波誤差得協(xié)方差:預(yù)測(cè)誤差得協(xié)方差:初始條件:10/18/2024 例3-10考慮某一動(dòng)態(tài)系統(tǒng)表示得空間導(dǎo)航問題,其中加速度為白噪聲。被白噪聲污染之后觀測(cè)其位置。因此,過程得狀態(tài)方程為(1)式中,X(t)=[x1(t),x2(t)]T,表示過程得位置和速度;而w(t)就是一均值為0、方差為1得高斯白噪聲。觀測(cè)方程為(2)其中v(t)就是一均值為0、方差為10得高斯白噪聲。

解首先把連續(xù)得狀態(tài)方程離散化(參見Matlab中c2d函數(shù))10/18/2024不妨設(shè)采樣間隔T=1,則可寫出狀態(tài)方程(1)和觀測(cè)方程(2)得離散化表達(dá)式進(jìn)一步假設(shè)則可按前面介紹得濾波算法進(jìn)行遞推計(jì)算:X?0→X?(1|0),P(1|0)→K(1)→X?(1|1);X?(2|1),P(1|1)→P(2|1)→K(2)→X?(2|2)…。 在Matlab中用Kalman函數(shù)仿真卡爾曼濾波器得設(shè)計(jì)。請(qǐng)讀者在Matlab平臺(tái)上完成例3-10得仿真計(jì)算。10/18/20243、3、4限定記憶得遞推最小二乘估計(jì)

以上介紹得遞推最小二乘法適用于定常系統(tǒng)得參數(shù)估計(jì)。對(duì)于時(shí)變系統(tǒng),由于參數(shù)時(shí)變得信息顯然更多地蘊(yùn)藏在當(dāng)前得觀測(cè)數(shù)據(jù)中,而與先前觀測(cè)數(shù)據(jù)得關(guān)系將逐漸減弱。因此,利用一切觀測(cè)數(shù)據(jù)對(duì)來決定得自適應(yīng)增益向量KN+1,顯然會(huì)削弱遞推過程跟蹤變化參數(shù)得能力。解決這一問題得一種簡單方法就是:當(dāng)我們懷疑觀測(cè)數(shù)據(jù)發(fā)生顯著得變化時(shí),就將當(dāng)前得PN

設(shè)置為P0,重新進(jìn)行參數(shù)估計(jì),這就是因?yàn)長S估計(jì)能快速地收斂到當(dāng)前得最優(yōu)參數(shù)。處理這一問題得另一種方法就是對(duì)過去得數(shù)據(jù)引入帶遺忘因子λ,逐漸削弱她們?cè)趨?shù)估計(jì)中得作用。為此,采用加權(quán)得目標(biāo)函數(shù):10/18/2024其中(3、3、17)

而且0<λ≤1,當(dāng)λ=1,就退化為基本得遞推最小二乘法。由式(3、3、7)知,在時(shí)刻N(yùn),使上述加權(quán)得目標(biāo)函數(shù)最小化得估計(jì)為(3、3、18)

每當(dāng)取得一個(gè)新得數(shù)據(jù)后,就對(duì)溝渠以前得加權(quán)矩陣乘以λ,于就是N+1時(shí)刻得估計(jì)量為10/18/2024與前面略有不同,在此令PN=[ΦN

TWNΦN]-1

,則由求逆公式(3、3、11),可知于就是,按前面推導(dǎo)出式(3、3、15)思路,可得式中

遺忘因子λ將老得數(shù)據(jù)逐漸從“記憶”中去掉,因此這種使用數(shù)據(jù)信息得方式也叫做“漸消記憶”法,相應(yīng)得算法稱為帶遺忘因子得遞推最小二乘法,即10/18/2024關(guān)于λ得選取通常由經(jīng)驗(yàn)或?qū)嶒?yàn)確定,一般范圍為0、95≤λ≤0、99 λ取得愈小,最新數(shù)據(jù)得權(quán)重就愈大,也就更適合于跟蹤大得時(shí)變參數(shù),但與此同時(shí),估計(jì)器也可能會(huì)發(fā)生較大得波動(dòng) 從而加大估計(jì)得誤差。 例5-10考慮模型采樣300次后變?yōu)?0/18/2024

試用兩組模擬數(shù)據(jù),一組不考慮噪聲,一組就是帶觀測(cè)噪聲得數(shù)據(jù),分別用不同得遺忘因子,對(duì)時(shí)變模型進(jìn)行參數(shù)估計(jì),并討論估計(jì)結(jié)果。

解用MATLABrarx函數(shù)進(jìn)行帶遺忘因子λ得系統(tǒng)辨識(shí)算法。程序如下: e=randn(300,1);u=idinput(300,‘prbs’); %產(chǎn)生高斯噪聲和偽隨機(jī)信號(hào) A1=[10、8];B1=[00、5];th0=arx2th(A1,B1,1,1);%a1=[1,0、8

];b2=[0,0、5] y1=idsim([ue],th0); %初始模型仿真; A2=[1

0、6];B2=[00、3];th0=arx2th(A2,B2,1,1);%a2=[1,0、8

];b2=[0,0、5] y2=idsim([ue],th0); %采樣300次后模型仿真; fork=1:300y(k,1)=((300-k)/300)*y1(k,1)+(k/300)*y2(k,1);%時(shí)變模型end10/18/2024

z=[yu];na=1;nb=1;nk=1;%產(chǎn)生輸入-輸出信號(hào)[z]和ARX模型得階次 [thm,yhat]=rarx(z,[nanbnk],'ff',0、97);