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

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

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

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

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

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

假定數(shù)據(jù){xk}就是獨立觀測得,容易驗證

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

(2)有效性 如果1

和2

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

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

稱為有效估計量。

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

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

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

和θ2

就是未知參數(shù)。

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

和2

分別就是θ1

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

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

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

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

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

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

(3、1、3)

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

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

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

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

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

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

)就是非負得,因此,?RM

MSE

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

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

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

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

獨立;(ⅱ)e與x

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

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

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

(y)。

解如果將x

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

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

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

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

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

(y)

(ii)當(dāng)e

與x

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

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

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

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

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

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

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

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

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

就是無偏得,而方差得ML估計

ML就是有偏得。但若將ML

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

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

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

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

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

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

信號和噪聲得協(xié)方差矩陣可分別表示為(3、2、7)10/18/2024于就是,觀測樣本得似然函數(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ù)據(jù)得長度),W就是接收過程得帶寬。 容易驗證,行列式|Cw

+Cs

|與時延τ無關(guān)。于就是,ML估計就就是選擇τ,使lnp(X|τ)最大,也即使式(3、2、8)得指數(shù)函數(shù)(3、2、9)10/18/2024最大。下面,我們從式(3、2、9)出發(fā),推導(dǎo)時延參數(shù)τ得最大似然估計得等效形式。為此,首先引進下列求逆公式(3、2、10)式中,A為n×n非奇異矩陣;g為n×1列向量。證明留給請讀者課外練習(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)最大,等價于使下式(3、2、11)

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

(3、2、12)

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

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

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

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

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

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

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

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

(3、3、2)其中

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

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

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

(3、3、4)

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

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

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

(3、3、5)

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

=LS。

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

和θ1得LS估計量。

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

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

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

[p,s]=polyfit(x,y,1) %生成擬合一次多項式運行結(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è)有一個非線性模型得輸出為其中,x1,x2

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

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

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

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

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

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); %實際系統(tǒng)得ARX模型 e=randn(200,1);u=idinput(200,‘prbs’);%高斯噪聲和偽隨機信號 y=idsim([ue],th0);z=[yu]; %模型仿真;輸入-輸出信號[z] na=2;nb=3;nk=1 %ARX模型得階次 [thm,yhat]=rarx(z,[nanbnk],'ng',0、1);%根據(jù)[z]進行ARX模型參數(shù)辨識 plot(y,'-');grid %作圖,實際系統(tǒng)得輸出曲線 holdon plot(yhat,':') %作圖,辨識系統(tǒng)得輸出曲線

參數(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)指出,要求觀測數(shù)據(jù)容量N≥n+m+1就是為了保證ΦTΦ非奇異,降低過程噪聲序列{e(k)}得影響,從而提高參數(shù)估計得精度。不論{e(k)}就是何種形式得噪聲序列,式(3、3、5)總就是成立得。換言之,噪聲性質(zhì)僅影響LS估計得統(tǒng)計特性。 下面介紹LS估計得統(tǒng)計特性。如果觀測噪聲或建模誤差序列{e(k)}具有零均值和相同得方差,即則LS估計量LS就是無偏、有效和相容得,并具估計誤差得協(xié)方差為σ

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

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

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

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

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

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

(3、3、11)

下面介紹最小二乘遞推算法。為了簡化符號,以下推導(dǎo)均用代替LS。 設(shè)ΦN就是時刻N為止得觀測數(shù)據(jù),N+1時刻θ得LS估計為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)

上式表明,新得估計量N+1等于前一時刻得估計量N

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

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

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

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

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

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

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

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

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

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

解首先把連續(xù)得狀態(tài)方程離散化(參見Matlab中c2d函數(shù))10/18/2024不妨設(shè)采樣間隔T=1,則可寫出狀態(tài)方程(1)和觀測方程(2)得離散化表達式進一步假設(shè)則可按前面介紹得濾波算法進行遞推計算: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è)計。請讀者在Matlab平臺上完成例3-10得仿真計算。10/18/20243、3、4限定記憶得遞推最小二乘估計

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

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

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

每當(dāng)取得一個新得數(shù)據(jù)后,就對溝渠以前得加權(quán)矩陣乘以λ,于就是N+1時刻得估計量為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)驗或?qū)嶒灤_定,一般范圍為0、95≤λ≤0、99 λ取得愈小,最新數(shù)據(jù)得權(quán)重就愈大,也就更適合于跟蹤大得時變參數(shù),但與此同時,估計器也可能會發(fā)生較大得波動 從而加大估計得誤差。 例5-10考慮模型采樣300次后變?yōu)?0/18/2024

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

解用MATLABrarx函數(shù)進行帶遺忘因子λ得系統(tǒng)辨識算法。程序如下: e=randn(300,1);u=idinput(300,‘prbs’); %產(chǎn)生高斯噪聲和偽隨機信號 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);%時變模型end10/18/2024

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