Lorenz系統(tǒng)參數(shù)估計(jì)方法研究

章翠蓮，李維德*，朱高峰

(1.蘭州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 甘肅 蘭州 730000; 2.蘭州大學(xué)資源環(huán)境學(xué)院, 甘肅 蘭州 730000)

Lorenz系統(tǒng)參數(shù)估計(jì)方法研究

章翠蓮1，李維德1*，朱高峰2

(1.蘭州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 甘肅 蘭州 730000; 2.蘭州大學(xué)資源環(huán)境學(xué)院, 甘肅 蘭州 730000)

馬爾科夫鏈蒙特卡洛方法(簡(jiǎn)稱MCMC)可用于復(fù)雜系統(tǒng)的不確定性估計(jì)和參數(shù)估計(jì).基于貝葉斯理論，運(yùn)用幾種改進(jìn)的MCMC方法：自適應(yīng)Metropolis 算法(AM)、延遲拒絕自適應(yīng)Metropolis算法(DRAM)、差分進(jìn)化馬爾科夫鏈算法(DE-MC)對(duì)具有復(fù)雜的動(dòng)態(tài)性質(zhì)的Lorenz 混沌系統(tǒng)未知參數(shù)進(jìn)行了探討性的估計(jì).根據(jù)未知參數(shù)的后驗(yàn)概率密度似然函數(shù)，利用MATLAB 仿真，選取樣本點(diǎn)出現(xiàn)頻率高的區(qū)間作為目標(biāo)分布區(qū)域，并通過(guò)縮小先驗(yàn)分布范圍來(lái)計(jì)算參數(shù)的估計(jì)值.分析比較這三個(gè)算法模擬的結(jié)果，得出如下結(jié)論：合適的目標(biāo)分布區(qū)域在參數(shù)估計(jì)中很關(guān)鍵；在搜索樣本點(diǎn)上，DRAM 算法的遍歷性最高；DE-MC 算法可使Lorenz 系統(tǒng)獲得較為精確的參數(shù)向量，更適用于復(fù)雜系統(tǒng)未知參數(shù)的估計(jì).

AM；DRAM；差分進(jìn)化馬爾科夫鏈；Lorenz 混沌系統(tǒng)；參數(shù)估計(jì)

1963年，洛倫茲提出了第一個(gè)表現(xiàn)奇異吸引子的動(dòng)力學(xué)系統(tǒng)，即Lorenz系統(tǒng)[1].該系統(tǒng)表現(xiàn)出非線性動(dòng)力學(xué)系統(tǒng)的復(fù)雜形式——具有混沌動(dòng)態(tài)性質(zhì).近年來(lái)，混沌系統(tǒng)在保密通信，人工智能以及生命科學(xué)等領(lǐng)域得到應(yīng)用和發(fā)展.但由于混沌系統(tǒng)的復(fù)雜性，參數(shù)的不可觀測(cè)性，或者保密性，系統(tǒng)的參數(shù)是未知的.因此系統(tǒng)的參數(shù)估計(jì)在混沌控制以及同步領(lǐng)域上的作用是不可估量的[2].

通過(guò)觀測(cè)數(shù)據(jù)求模型的未知參數(shù)，是建模分析中的重要步驟，但由于觀測(cè)過(guò)程中存在系統(tǒng)誤差以及觀測(cè)信息的偶然誤差，不一定可求得合適的參數(shù)解，或者得出的參數(shù)解不唯一，因此本文采用基于貝葉斯理論下改進(jìn)的MCMC方法對(duì)Lorenz系統(tǒng)未知參數(shù)估計(jì)進(jìn)行探討性的研究.

MCMC的Metropolis-Hasting(M-H)算法，其思想是構(gòu)造一個(gè)以目標(biāo)分布π(x)為不變分布的馬爾科夫鏈,通過(guò)一個(gè)概率密度函數(shù)q(x,y) (稱為建議函數(shù))來(lái)實(shí)現(xiàn)這一分布[3,4]. 該方法是否具有收斂性取決于是否選取合適的建議分布函數(shù).假設(shè)當(dāng)前狀態(tài)為x，則從建議函數(shù)q(x,.)產(chǎn)生一個(gè)候選點(diǎn)y，那么接受概率:

(1)

若α>u(u是[0,1]均勻分布的隨機(jī)數(shù))，接受候選點(diǎn)y作為新的采樣點(diǎn)，否則復(fù)制原來(lái)的x作為新的采樣點(diǎn).

當(dāng)建議分布為對(duì)稱分布時(shí)，即q(x,y)=q(y,x).此時(shí)接受概率為：

(2)

近年來(lái)，許多研究人員在MCMC算法的運(yùn)行效率和效果做出許多努力從而改進(jìn)了方法.HeikkiHaario等提出了自適應(yīng)Metropolis算法，簡(jiǎn)稱AM.這個(gè)算法利用了MarkovChain的歷史信息自動(dòng)調(diào)節(jié)提議函數(shù)的協(xié)方差矩陣，與隨機(jī)游走算法相比，提高了算法運(yùn)行效率[5,6].文獻(xiàn)[2] 應(yīng)用了該算法對(duì)Lorenz系統(tǒng)的參數(shù)進(jìn)行了估計(jì).在AM算法提出后，HeikkiHaario等又提出了自適應(yīng)Metropolis采樣和延遲拒絕相結(jié)合的算法：延遲拒絕自適應(yīng)Metropolis算法，簡(jiǎn)稱DRAM.這個(gè)算法提高了馬爾科夫鏈的遍歷性和采樣點(diǎn)的接受率[7].TerBraak則通過(guò)遺傳算法思想，提出了一種差分進(jìn)化馬爾科夫鏈算法，簡(jiǎn)稱DE-MC.該算法多條鏈同時(shí)運(yùn)行，其解決了MCMC的一個(gè)重要問(wèn)題，即為跳躍分布提供了合適的范圍和方向[8].基于貝葉斯理論，本文主要將這三種改進(jìn)的MCMC算法運(yùn)用到對(duì)Lorenz系統(tǒng)的未知參數(shù)估計(jì)，對(duì)這三種方法進(jìn)行探討性的討論，得出解決Lorenz系統(tǒng)未知參數(shù)估計(jì)的較優(yōu)方法.

1 理論與算法

1.1 貝葉斯理論

貝葉斯方法提供了一個(gè)涵蓋模型不確定性的框架，在這個(gè)框架下，具有概率分布的貝葉斯方法的關(guān)鍵的特點(diǎn)在于描述了參數(shù)和模型的不確定性[9].

基于貝葉斯理論，我們認(rèn)為未知參數(shù)向量(對(duì)Lorenz系統(tǒng)而言，其有3個(gè)未知參數(shù)，因此d=3),在對(duì)X進(jìn)行初步推斷時(shí)，根據(jù)相關(guān)未知參數(shù)的信息 ，我們認(rèn)為其先驗(yàn)分布為p(θ).在對(duì)系統(tǒng)進(jìn)行觀測(cè)采樣后，得到觀測(cè)信息(包含N個(gè)樣本點(diǎn)).因此，由貝葉斯估計(jì)，系統(tǒng)參數(shù)的后驗(yàn)分布，觀測(cè)信息和先驗(yàn)分布之間的關(guān)系[2]：

(3)

考慮到Z是已知的信息，上式可以寫(xiě)成：

p(θ|Z)∝p(Z|θ)P(θ)

(4)

我們知道，貝葉斯推斷的核心在于使用似然函數(shù)去分析參數(shù)的不確定性.因此對(duì)p(Z|), 在實(shí)際中一般用似然函數(shù)表示，其值越大表明擬合程度越好[10].

1.2AM算法

AM是基于Metropolis算法中的隨機(jī)游走算法(RWM)[11]改進(jìn)而來(lái)的,該算法關(guān)鍵在于，先構(gòu)造具有高斯分布的建議函數(shù)，利用馬爾科夫鏈先前全部樣本信息計(jì)算建議分布的協(xié)方差矩陣，自適應(yīng)地向目標(biāo)函數(shù)逼近[12]. 利用該點(diǎn)，通過(guò)后驗(yàn)分布獲得的信息使得建議分布得到更新.在第i步中，Haario等提出了以當(dāng)前點(diǎn)作為均值和具有協(xié)方差矩陣Ci的多元正態(tài)分布.該協(xié)方差Ci的定義為：

(5)

其中，sd是縮放因子，a1>0 是一個(gè)很小的數(shù)，它們分別取為2.382/d與10-6.在理論上要保證Ci非奇異，Id是d維單位矩陣.C0是初始協(xié)方差矩陣.Cov(θ0,...,θi-1)表示θ0,...,θi-1的協(xié)方差矩陣. 若初始協(xié)方差過(guò)大也會(huì)使AM算法自適應(yīng)性變差；反之，初始協(xié)方差過(guò)小降低AM算法的遍歷性. 為非自適應(yīng)長(zhǎng)度，其值越大自適應(yīng)性越慢[12].

AM算法流程如下：

a)設(shè)定i=1，對(duì)所求的參數(shù)變量初始化；

b)利用(5)式計(jì)算協(xié)方差矩陣Ci；

c)產(chǎn)生候選點(diǎn)θ*～N(θi-1,Ci)；

d)根據(jù)(1)式，計(jì)算接受概率α=min{1,p(θ*|Z)/p(θi-1|Z)}；

e)產(chǎn)生滿足均勻分布的隨機(jī)數(shù)u～U[0,1]；

f)若α>u,則θi=θ*，否則θi=θi-1；

g)重復(fù)b)～f)，直到完成產(chǎn)生指定的迭代數(shù)量的樣本.

1.3DRAM算法

DRAM是基于延遲拒絕算法[13]和AM算法相結(jié)合而得來(lái)的算法[7].

DR算法是具有不同提議函數(shù)和轉(zhuǎn)移核的MCMC方法.Peskun和Tierney已證明出有限維的狀態(tài)空間和普遍的狀態(tài)空間里DR算法通過(guò)降低停留在當(dāng)前候選點(diǎn)的概率提高M(jìn)CMC方法的運(yùn)行效率[14,15].其主要思想在于：在進(jìn)行M-H抽樣時(shí)，若當(dāng)前的候選點(diǎn)被拒絕，不是保持當(dāng)前狀態(tài)不變，而是從當(dāng)前位置利用改變的提議函數(shù)再次移動(dòng)，通過(guò)對(duì)第2次移動(dòng)的候選點(diǎn)的接受概率計(jì)算來(lái)保持馬爾科夫鏈的可逆性，其本質(zhì)上是通過(guò)對(duì)提議函數(shù)局部自適應(yīng)調(diào)整來(lái)提高算法的效率[12].

在算法進(jìn)程中，假定當(dāng)前位置為x，用π(·)表示該點(diǎn)的目標(biāo)分布函數(shù)，根據(jù)多元正態(tài)分布函數(shù)q0(x,·)產(chǎn)生新的候選點(diǎn)y，接受率(類似式(1))為：

(6)

若在點(diǎn)x被拒絕，則根據(jù)當(dāng)前位置和被拒的候選點(diǎn)再次移動(dòng)，從新的建議函數(shù)q1(x,y,·)產(chǎn)生新的候選點(diǎn)z,相應(yīng)的接受概率為[16]：

(7)

若是第二級(jí)的候選點(diǎn)仍被拒絕，則可在被拒的候選點(diǎn)z上再進(jìn)行第三級(jí)移動(dòng).

DRAM與AM算法的流程大致相同，不同于AM算法的是在第f)步.將第f)步改動(dòng):若α>u, 則θi=θ*，再返回到第b)步繼續(xù)迭代；若α

1.4DE-MC算法

差分進(jìn)化馬爾科夫鏈?zhǔn)窃诰哂蠱CMC模擬的實(shí)參數(shù)空間上利用Metropolis準(zhǔn)則與差分進(jìn)化算法[17]結(jié)合而得來(lái)的MCMC算法[8].該算法克服了MCMC方法中建議分布需要選取合適的范圍和方向以及相應(yīng)的計(jì)算效率問(wèn)題.該算法有兩個(gè)參數(shù)：縮放因子γ和平行鏈的數(shù)目N可根據(jù)實(shí)際情況定義[18].

DE-MC算法的基本過(guò)程[8]：

θ*=θi+γ(θa-θb)+ε

(8)

其中，θ*是建議選取的參數(shù)集，θi是當(dāng)前的參數(shù)集，θa,θb是從θi不含外的參數(shù)集隨機(jī)選取出來(lái)的，γ為縮放因子.ε～N(0,b)d,在這里b很小.建議的參數(shù)集接受或者拒絕取決于Metropolis比率,具體可參考(1) 式[16].

上面所列舉的三種MCMC方法，AM和DRAM算法都具有自適應(yīng)性，DE-MC是遺傳算法和傳統(tǒng)的Metropolis方法相結(jié)合的算法.每一種算法各有其相關(guān)的優(yōu)勢(shì)，因此下面將這三種算法運(yùn)用到Lorenz模型進(jìn)行未知參數(shù)的估計(jì)，進(jìn)行探討性的討論.

2 Lorenz混沌系統(tǒng)參數(shù)估計(jì)分析

Lorenz系統(tǒng)可由如下的常微分方程組表示：

(9)

我們知道當(dāng)系統(tǒng)(9)中a,b,r分別為10,28,8/3時(shí)，系統(tǒng)表現(xiàn)為混沌現(xiàn)象.Lorenz描述的非線性耗散系統(tǒng)具有初始條件敏感性.在用數(shù)值解法求解非線性微分問(wèn)題時(shí)，數(shù)值積分方法和編程現(xiàn)實(shí)中的誤差累積導(dǎo)致模擬的軌道與原來(lái)的軌道偏離.即不同的初始條件，不同的算法和編程細(xì)節(jié)，使得計(jì)算軌道是不可重復(fù)的.因此，系統(tǒng)的長(zhǎng)期行為是不可預(yù)測(cè)的.本文選取的時(shí)間長(zhǎng)度為10s.首先對(duì)該系統(tǒng)進(jìn)行數(shù)據(jù)采樣：先讓Lorenz系統(tǒng)進(jìn)行演化,以混沌系統(tǒng)中某一點(diǎn)(1.500，2.230，-1.270)作為初始點(diǎn)，并記為：(x1(0),x2(0),x3(0)).設(shè)置時(shí)間步長(zhǎng)h為0.01s,采用四階定步長(zhǎng)Runge-Kutta算法來(lái)求解常微分方程，得出離散的數(shù)值解，我們得到Lorenz系統(tǒng)在10s內(nèi)的1000組離散數(shù)據(jù)： (x1(ti),x2(ti),x3(ti))(i=0,1,...,1000).

它們可作為觀測(cè)數(shù)據(jù)集Z, 由于這些數(shù)據(jù)在操作過(guò)程中必有觀測(cè)誤差，因此可多模擬幾次，取每個(gè)步長(zhǎng)狀態(tài)的平均值作為觀測(cè)數(shù)據(jù)集來(lái)減少誤差的干擾.

現(xiàn)在已知觀測(cè)數(shù)據(jù)集，分別采用AM,DRAM和DE-MC估計(jì)出Lorenz系統(tǒng)的未知參數(shù)集θ=(a,b,r).由貝葉斯公式求出參數(shù)的聯(lián)合后驗(yàn)概率密度函數(shù).令模型的初始狀態(tài)量為：(x1(0),x2(0),x3(0))，每隔10個(gè)步長(zhǎng)共取10個(gè)樣本點(diǎn)，即抽取0h，10h，...，90h時(shí)刻的狀態(tài)量(x1(ti),x2(ti),x3(ti))(i=1,2,...,10)作為這次的觀測(cè)數(shù)據(jù)集Z.

根據(jù)貝葉斯公式(3),未知參數(shù)的后驗(yàn)分布：

p(a,b,r|Z)=k·p(Z|a,b,r)p(a,b,r)

(10)

其中，k是一個(gè)大于0的常數(shù).

在貝葉斯推斷中，首先設(shè)定未知參數(shù)的先驗(yàn)分布p(a,b,r).假設(shè)這三個(gè)參數(shù)均滿足獨(dú)立均勻分布，用U(a),U(b),U(r)分別表示三個(gè)參數(shù)的均勻分布，則

p(a,b,r)=U(a)U(b)U(r)

(11)

相比較其他分布，均勻分布是最簡(jiǎn)單的分布，其限定了參數(shù)變化的界限，先驗(yàn)分布的界限很小時(shí)，目標(biāo)區(qū)域則相對(duì)變小，這樣有利于提高參數(shù)的精確度.通常情況下熟知?dú)v史經(jīng)驗(yàn)更能得到精確的先驗(yàn)分布.

對(duì)于式(10)中的p(Z|a,b,r),前面提到觀測(cè)集是經(jīng)過(guò)多次測(cè)量取平均值而得到的.根據(jù)普適似然不確定性分析方法GLUE的基本原理,用似然函數(shù)來(lái)比較模擬值和實(shí)測(cè)值的擬合程度，本文選用基于殘差和的似然函數(shù)[19,20]：

(12)

本文取n=10,d=3,M=1.對(duì)向量進(jìn)行更新后，通過(guò)用四階定步長(zhǎng)Runge-Kutta算法解Lorenz系統(tǒng)微分方程，時(shí)間和步長(zhǎng)不變，每隔10h取一個(gè)數(shù)據(jù)，則得到10個(gè)數(shù)據(jù).第i個(gè)數(shù)據(jù)記為:(y1(ti),y2(ti),y3(ti))(i=1,2,...,10).

由(10)-(12)可得，該系統(tǒng)未知參數(shù)的后驗(yàn)概率密度似然函數(shù)為：

(13)

在模擬之前，本文先假設(shè)三個(gè)參數(shù)的目標(biāo)分布區(qū)間均為[0,30],即a,b,r的聯(lián)合先驗(yàn)分布均為：

(14)

后驗(yàn)概率密度似然函數(shù)確定后，則可以用AM,DRAM,DE-MC算法對(duì)這三個(gè)未知參數(shù)進(jìn)行隨機(jī)抽樣.對(duì)于AM和DRAM算法,非自適應(yīng)段長(zhǎng)度=200.DRAM算法中，λ=0.1，m=3.

3 仿真結(jié)果及分析

考慮到三個(gè)參數(shù)的目標(biāo)分布都較大，迭代次數(shù)先取為2000，參數(shù)初始值則取相對(duì)應(yīng)目標(biāo)分布中的隨機(jī)值. 在每個(gè)算法進(jìn)行抽樣時(shí)，每次都對(duì)未知參數(shù)向量進(jìn)行更新，通過(guò)計(jì)算該隨機(jī)參數(shù)向量的似然后驗(yàn)密度函數(shù)值與上一步后驗(yàn)密度函數(shù)值的比值，決定是否將該向量作為下一步的值并進(jìn)行更新.

在各參數(shù)的目標(biāo)分布區(qū)間皆為[0,30]時(shí)，各個(gè)算法模擬出來(lái)各參數(shù)的馬爾科夫鏈還未收斂(圖1)，仿真的結(jié)果不理想.根據(jù)圖1中基于殘差平方和似然函數(shù)值較大的部分對(duì)應(yīng)的參數(shù)向量，即a=9.9529,b=28.2056,r=2.7923時(shí)，式(12)中的較大，說(shuō)明我們估計(jì)的參數(shù)向量在這個(gè)向量的附近.但由于MCMC算法具有隨機(jī)性，三種算法模擬出來(lái)三個(gè)參數(shù)的直方圖都不理想(圖2)， 縮小各個(gè)參數(shù)的目標(biāo)分布區(qū)域再進(jìn)行Matlab模擬可期望得出更精確的參數(shù)向量.重新設(shè)置參數(shù)的目標(biāo)分布：a～U[9.9529-2,9.9529+2],b～U[28.2056-2,30],r～U[2.7923-2，2.7923+2]，其他條件不變.

圖1 目標(biāo)分布區(qū)間為[0,30]下各算法仿真結(jié)果

注：T表示迭代次數(shù)，Q表示目標(biāo)分布區(qū)域，L表示基于殘差和的似然函數(shù)值；(1)-(3)分別是用AM模擬得的參數(shù)a,b,r的馬爾科夫鏈；(5)-(7)分別是用$DRAM$模擬得的參數(shù)a,b,r的馬爾科夫鏈；(9)-(11)分別是用DE-MC模擬得參數(shù)a,b,r的馬爾科夫鏈；(4),(8),(12)分別是用AM,DRAM,DE-MC模擬得到的基于殘差和的似然函數(shù)值演化過(guò)程.

圖2 目標(biāo)分布區(qū)間為[0,30]下各算法仿真的參數(shù)值方圖

注：P表示頻率；(1)-(3)分別是用AM模擬出來(lái)參數(shù)a,b,r的直方圖；(4)-(6)分別是用DRAM模擬出來(lái)參數(shù)a,b,r的直方圖；(7)-(9)分別是用DE-MC模擬出來(lái)參數(shù)a,b,r的直方圖.

圖3 第一次調(diào)整目標(biāo)分布后各算法仿真結(jié)果

在新的參數(shù)目標(biāo)分布區(qū)域下，由三個(gè)算法模擬出各參數(shù)的馬爾科夫鏈以及計(jì)算出式(12)的值的演化過(guò)程(圖3)，我們可以看出DE-MC算法模擬出來(lái)的參數(shù)r的馬爾科夫鏈已收斂，其他兩個(gè)參數(shù)在小范圍里波動(dòng)，其基于殘差的平方和的似然函數(shù)值大體比其他兩個(gè)算法的大.在迭代次數(shù)1700-1800范圍內(nèi)，式(12)的值最大，此時(shí)參數(shù)向量為(9.9875,28.1047,2.6508)，AM和DRAM算法模擬各參數(shù)的馬爾科夫鏈皆無(wú)收斂的趨勢(shì).從遍歷性來(lái)看，DRAM的遍歷性最好(候選點(diǎn)接受率約為0.5)，AM算法次之(候選點(diǎn)接受率約為0.35)，DE-MC的遍歷性最差.結(jié)合三個(gè)算法模擬出各個(gè)參數(shù)的直方圖(圖4)，無(wú)法從AM算法模擬出樣本點(diǎn)的直方圖減小目標(biāo)分布區(qū)域；DRAM算法模擬出來(lái)的各參數(shù)樣本點(diǎn)直方圖類似正態(tài)分布，參數(shù)r的最為明顯；從DE-MC算法模擬出的各參數(shù)樣本點(diǎn)直方圖可以看出樣本點(diǎn)皆分布集中，因此可根據(jù)其縮小目標(biāo)分布區(qū)域.考慮DE-MC算法中的式(12)似然函數(shù)較大值所對(duì)應(yīng)的參數(shù)向量，可將參數(shù)向量的目標(biāo)分布設(shè)置為：(9.9875±0.5,28.1047±0.5,2.6508 ±0.2).

圖4 第一次目標(biāo)分布調(diào)整后各算法仿真的參數(shù)值方圖

利用調(diào)整后的目標(biāo)分布區(qū)域得出三個(gè)算法模擬出來(lái)的各參數(shù)的馬爾科夫鏈以及由式(12)得到的似然函數(shù)值演化圖(圖5).從圖5可看出，DRAM算法與其他兩個(gè)算法相比，其遍歷性最高，樣本點(diǎn)的接受概率約為0.54,AM算法樣本點(diǎn)的接受概率越為0.14，但是這兩個(gè)算法最后未收斂.從后驗(yàn)概率密度似然函數(shù)來(lái)看，DE-MC算法收斂，并且其似然函數(shù)值最大.選取似然函數(shù)值最大的部分，即用DE-MC算法模擬出的馬爾科夫鏈?zhǔn)諗康哪遣糠?，其?duì)應(yīng)的參數(shù)向量約為：(9.9949, 28.0042, 2.6673), 即約為9.9949，約為28.0042，約為2.6673. 將其帶入(9)式，經(jīng)過(guò)數(shù)值模擬，由(12)式可得出 約為1006，即其與觀測(cè)集的殘差平方和約為0.001. 由此可見(jiàn)，DE-MC算法更適用于Lorenz混沌系統(tǒng)的未知參數(shù)估計(jì).

圖5 第二次目標(biāo)分布調(diào)整后各算法仿真結(jié)果

圖6 狀態(tài)變量估計(jì)值與真實(shí)值在15-25s內(nèi)的軌跡對(duì)比圖

Matlab中解一階微分方程組初值問(wèn)題常見(jiàn)的方法為ode45函數(shù)，該函數(shù)實(shí)現(xiàn)了變步長(zhǎng)四階五級(jí)的Runge-Kutta-Felhberg算法，采用變步長(zhǎng)算法解微分方程[21].將原來(lái)的參數(shù)向量和通過(guò)DE-MC算法求解出來(lái)的參數(shù)向量代入Lorenz模型的狀態(tài)方程(9)，利用MATLAB中ode45函數(shù)數(shù)值解法求解.當(dāng)時(shí)間長(zhǎng)度為10s，初始狀態(tài)相同，發(fā)現(xiàn)各狀態(tài)變量真實(shí)值和估計(jì)值在10s內(nèi)的振幅以及頻率相似，相位大體相同，軌跡幾乎一模一樣.其他條件不變，延長(zhǎng)時(shí)間長(zhǎng)度，各狀態(tài)變量真實(shí)值和估計(jì)值在20s后開(kāi)始不在一條軌道上，估計(jì)值漸漸脫離原來(lái)的軌道(圖6).因此，Lorenz系統(tǒng)的長(zhǎng)期行為是不可預(yù)測(cè)的.在選取觀測(cè)樣本數(shù)據(jù)時(shí)，系統(tǒng)不宜長(zhǎng)時(shí)間演化，時(shí)間一般不超過(guò)10s. 在適當(dāng)?shù)哪繕?biāo)分布區(qū)域內(nèi)，DE-MC算法適用于短時(shí)間演化的復(fù)雜系統(tǒng)的參數(shù)估計(jì)，其最后模擬出的參數(shù)向量是有意義的.

4 結(jié)論與討論

本文基于貝葉斯理論的框架，利用MCMC方法中的AM,DRAM和DE-MC算法對(duì)Lorenz模型未知參數(shù)向量進(jìn)行采樣，經(jīng)過(guò)兩次目標(biāo)分布區(qū)域的調(diào)整，分析比較各參數(shù)的馬爾科夫鏈，得到了系統(tǒng)的參數(shù)估計(jì)值.本文研究的主要結(jié)論總結(jié)如下：

(1)當(dāng)未知參數(shù)的目標(biāo)分布區(qū)域較大時(shí)，這三個(gè)算法模擬出各參數(shù)的馬爾科夫鏈皆不收斂，基于殘差和的似然函數(shù)值不高，效果不是很理想.參照各算法模擬出的似然函數(shù)值演化過(guò)程圖以及各參數(shù)樣本點(diǎn)的直方圖，選取樣本點(diǎn)出現(xiàn)頻率高的區(qū)間作為新的目標(biāo)分布區(qū)域，及時(shí)調(diào)整先驗(yàn)分布.參數(shù)的先驗(yàn)分布越準(zhǔn)確，即目標(biāo)分布區(qū)域越小，模擬出來(lái)的參數(shù)向量越精確.

(2)本文應(yīng)用了MCMC的三種算法對(duì)Lorenz混沌模型未知參數(shù)進(jìn)行了估計(jì).與文獻(xiàn)[2]仿真的結(jié)果不同的是，從運(yùn)行過(guò)程來(lái)看，在相同的迭代次數(shù)下，AM算法運(yùn)行的時(shí)間最少；由于DE-MC算法是多條鏈同時(shí)運(yùn)行，模擬過(guò)程最緩慢.從樣本點(diǎn)的搜索區(qū)域來(lái)看，DRAM算法搜索范圍大，遍歷性最好；而DE-MC算法的遍歷性較差.從模擬出來(lái)的結(jié)果來(lái)看，經(jīng)過(guò)兩次目標(biāo)分布區(qū)域的調(diào)整，在三種算法中，DE-MC算法模擬出的馬爾科夫鏈?zhǔn)諗浚⑶一跉埐詈偷乃迫缓瘮?shù)值較大，說(shuō)明模擬出來(lái)的參數(shù)向量與真實(shí)的參數(shù)向量相差甚小.

(3)從參數(shù)估計(jì)的效果看，在適當(dāng)?shù)哪繕?biāo)分布區(qū)域內(nèi)，對(duì)于短時(shí)間演化的Lorenz混沌系統(tǒng)，DE-MC算法能較快地模擬出該參數(shù)向量.對(duì)一般的非線性以及動(dòng)態(tài)的模型，DE-MC算法同樣能估計(jì)模型的參數(shù).

綜上所述，AM,DRAM和DE-MC算法皆是通過(guò)構(gòu)造馬爾科夫鏈進(jìn)行隨機(jī)模擬，在對(duì)Lorenz混沌系統(tǒng)進(jìn)行未知參數(shù)估計(jì)的過(guò)程中，實(shí)際就是在先驗(yàn)分布界定的目標(biāo)區(qū)域?qū)?shù)向量進(jìn)行隨機(jī)且較好抽取樣本點(diǎn)的過(guò)程，也就是構(gòu)造馬爾科夫鏈的過(guò)程.通過(guò)Matlab仿真模擬，DE-MC算法能較好地對(duì)Lorenz混沌系統(tǒng)未知參數(shù)進(jìn)行估計(jì)，因此該算法能更好地應(yīng)用在復(fù)雜系統(tǒng)的參數(shù)估計(jì)上.

[1]劉福才, 梁曉明,Hénon. 混沌系統(tǒng)的廣義預(yù)測(cè)控制與同步快速算法[J].物理學(xué)報(bào),2005, (10):4584-4589.

[2]曹小群,宋君強(qiáng).基于MCMC方法的Lorenz混沌系統(tǒng)的參數(shù)估計(jì)[J].國(guó)防科技大學(xué)學(xué)報(bào),2010,(2):68-72.

[3]ChristropheA,NandoDF.AnintroductiontoMCMCformachinelearning[J].MachineLearning,2003,50,5-43.

[4]陳平，徐若曦.Metropolis-Hastings自適應(yīng)算法及其應(yīng)用[J].系統(tǒng)工程理論與實(shí)踐,2008，(1):100-108.

[5]HaarioH,SaksmanE,TamminenJ.AnadaptiveMetropolisalgoritnm[J].Bernoulli,2001,(2):223-242.

[6]SaksmanE,ViholaM.OntheergodicityoftheadaptiveMetropolisalgorithmonunboundeddomain[J].TheAnnalsofAppliedProbability,2010,(6):2178-2203.

[7]HaarioH,LaineM,MiraA,etal.DRAM:EfficientadaptiveMCMC[J].Statistics&Computing, 2006,(4):339-354.

[8]BraakCJ.AMarkovChainMonteCarloversionofthegeneticalgorithmdifferentialevolution:EasyBayesiancomputingforrealparameterspaces[J].Statistics&Computing, 2006, (16):239-249.

[9]LombardiMJ.Bayesianinferenceforstableditributions:ArandomwalkMCMCapproach[J].ComputationalStatistics&DataAnalysis,2007,(5):2688-2700.

[10]MarshallL,NottD,SharmaA.AcomparativestudyofMarkovchainMonteCarlomethodsforconceptualrainfall-runoffmodeling[J].WaterResourcesResearch, 2004, (40):183-188.

[11]KuczeraG,ParentE.MonteCarloassessmentofparameteruncertaintyinconceptualcatchmentmodels:TheMetropolisalgorithm[J].JournalofHydrology, 1998, (1-4):69-85.

[12]郝燕玲,單志明.基于DRAM算法的α穩(wěn)定分布參數(shù)估計(jì)[J].華中科技大學(xué)學(xué)報(bào)，2011,(10):73-78.

[13]TierneyL,MiraA.SomeadaptiveMonteCarlomethodsforBayesianinference[J].StatisticsinMedicine, 1998, (17-18):2507-2515.

[14]PeskunPH.OptimumMonte-CarlosamplingusingMarkovchains[J].BiometrikaTrust, 1973,(3):607-612.

[15]TierneyL.AnoteonMetropolis-Hastingskernelsforgeneralstatespaces[J].AnnalsofAppliedProbability,1998, (1):1-9.

[16]SmithTJ,MarshallLA.Bayesianmethodsinhydrologicmodeling:AstudyofrecentadvancementsinMarkovchainMonteCarlotechniques[J].WaterResourcesResearch, 2008, (12):67-76.

[17]StornR,PriceK.Differentialevolution-Asimpleandefficientheuristicforglobaloptimizationovercontinuousspaces[J].JournalofGlobalOptimization, 1997,(4):341-359.

[18]BraakCJFT,VrugtJA.DifferentialevolutionMarkovchainwithsnookerupdaterandfewerchains[J].Statistics&Computing, 2008,(4):435-446.

[19]BevenKJBinleyAM.Thefutureofdistributedmodels:Modelscalibrationanduncertaintyprediction[J].HydrologicProcess,1992,(2):279-298.

[20]FreerJK,BevenKJ.Bayesianestimationofuncertaintyinrunoffpredictionandthevalueofdata:AnapplicationoftheGLUEapproach[J].WaterResourcesResearch,1996,(7):2161-2173.

[21]薛定宇,陳陽(yáng)泉.高等應(yīng)用數(shù)學(xué)問(wèn)題的MATLAB求解(第3版)[M].北京：清華大學(xué)出版社，2013.

(責(zé)任編校：晴川)

Study on Methods of Parameters Estimation of Lorenz System

ZHANG Cuilian1,LI Weide1*, ZHU Gaofeng2

(1.School of Mathematics and Statistics, Lanzhou University, Lanzhou Gansu 730000, China; 2. School of Resources and Environment, Lanzhou University, Lanzhou Gansu 730000, China)

Markov Chain Monte Carlo method (in short as MCMC) is useful to the uncertainty and parameter estimation for complicated system. Based on the Bayesian theory, the paper proposes an explorative study by utilizing several MCMC methods, such as Adaptive Metropolis (AM), Delay Rejection Adaptive Metropolis (DRAM) and Differential Evolution Markov Chain (DE-MC) algorithm, to estimate the parameters of Lorenz system which has complicated dynamical property. According to the posterior probability density likelihood function of unknown parameter, utilizing MATLAB simulation, selecting the region where sampling points appearing with high frequency as target distribution area, and narrowing the scope of the prior distribution, we calculate the estimation value of parameters. By analyzing and comparing the simulated results of the three algorithms mentioned above, we achieved these conclusions: suitable target distribution area is critical to estimate parameter; DRAM algorithm has high ergodicity in searching sampling points; DE-MC algorithm is more appropriate for complicated system with unknown parameter estimation, which makes Lorenz system obtain precise parameter vector.

AM; DRAM; differential evolution Markov Chain; Lorenz chaotic system; parameter estimation

2016-10-14

國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào)：41571016)資助項(xiàng)目.

章翠蓮(1991— )，女，廣西欽州人，蘭州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院碩士生.研究方向：應(yīng)用數(shù)學(xué).

O212.1

A

1008-4681(2017)02-0005-06

*通訊作者：李維德(1967— )，男，甘肅蘭州人，蘭州大學(xué)教學(xué)與統(tǒng)計(jì)學(xué)院教授，博士.研究方向：數(shù)字生態(tài)學(xué).