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

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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

1 理論與算法

1.1 貝葉斯理論

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

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

(3)

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

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

(4)

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

1.2AM算法

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

(5)

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

AM算法流程如下：

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

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

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

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

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

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

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

1.3DRAM算法

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

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

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

(6)

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

(7)

若是第二級的候選點仍被拒絕，則可在被拒的候選點z上再進行第三級移動.

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

1.4DE-MC算法

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

DE-MC算法的基本過程[8]：

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

(8)

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

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

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

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

(9)

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

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

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

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

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

(10)

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

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

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

(11)

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

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

(12)

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

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

(13)

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

(14)

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

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

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

在各參數(shù)的目標(biāo)分布區(qū)間皆為[0,30]時，各個算法模擬出來各參數(shù)的馬爾科夫鏈還未收斂(圖1)，仿真的結(jié)果不理想.根據(jù)圖1中基于殘差平方和似然函數(shù)值較大的部分對應(yīng)的參數(shù)向量，即a=9.9529,b=28.2056,r=2.7923時，式(12)中的較大，說明我們估計的參數(shù)向量在這個向量的附近.但由于MCMC算法具有隨機性，三種算法模擬出來三個參數(shù)的直方圖都不理想(圖2)， 縮小各個參數(shù)的目標(biāo)分布區(qū)域再進行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ù)值演化過程.

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

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

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

在新的參數(shù)目標(biāo)分布區(qū)域下，由三個算法模擬出各參數(shù)的馬爾科夫鏈以及計算出式(12)的值的演化過程(圖3)，我們可以看出DE-MC算法模擬出來的參數(shù)r的馬爾科夫鏈已收斂，其他兩個參數(shù)在小范圍里波動，其基于殘差的平方和的似然函數(shù)值大體比其他兩個算法的大.在迭代次數(shù)1700-1800范圍內(nèi)，式(12)的值最大，此時參數(shù)向量為(9.9875,28.1047,2.6508)，AM和DRAM算法模擬各參數(shù)的馬爾科夫鏈皆無收斂的趨勢.從遍歷性來看，DRAM的遍歷性最好(候選點接受率約為0.5)，AM算法次之(候選點接受率約為0.35)，DE-MC的遍歷性最差.結(jié)合三個算法模擬出各個參數(shù)的直方圖(圖4)，無法從AM算法模擬出樣本點的直方圖減小目標(biāo)分布區(qū)域；DRAM算法模擬出來的各參數(shù)樣本點直方圖類似正態(tài)分布，參數(shù)r的最為明顯；從DE-MC算法模擬出的各參數(shù)樣本點直方圖可以看出樣本點皆分布集中，因此可根據(jù)其縮小目標(biāo)分布區(qū)域.考慮DE-MC算法中的式(12)似然函數(shù)較大值所對應(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ū)域得出三個算法模擬出來的各參數(shù)的馬爾科夫鏈以及由式(12)得到的似然函數(shù)值演化圖(圖5).從圖5可看出，DRAM算法與其他兩個算法相比，其遍歷性最高，樣本點的接受概率約為0.54,AM算法樣本點的接受概率越為0.14，但是這兩個算法最后未收斂.從后驗概率密度似然函數(shù)來看，DE-MC算法收斂，并且其似然函數(shù)值最大.選取似然函數(shù)值最大的部分，即用DE-MC算法模擬出的馬爾科夫鏈?zhǔn)諗康哪遣糠?，其對?yīng)的參數(shù)向量約為：(9.9949, 28.0042, 2.6673), 即約為9.9949，約為28.0042，約為2.6673. 將其帶入(9)式，經(jīng)過數(shù)值模擬，由(12)式可得出 約為1006，即其與觀測集的殘差平方和約為0.001. 由此可見，DE-MC算法更適用于Lorenz混沌系統(tǒng)的未知參數(shù)估計.

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

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

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

4 結(jié)論與討論

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

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

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

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

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

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

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

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

[4]陳平，徐若曦.Metropolis-Hastings自適應(yīng)算法及其應(yīng)用[J].系統(tǒng)工程理論與實踐,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].華中科技大學(xué)學(xué)報，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)用數(shù)學(xué)問題的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

國家自然科學(xué)基金(批準(zhǔn)號：41571016)資助項目.

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

O212.1

A

1008-4681(2017)02-0005-06

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