李 揚(yáng)

（山西省水利水電科學(xué)研究院 山西太原 030002）

0 引言

水文頻率分析工作是水利工程規(guī)劃設(shè)計(jì)的基礎(chǔ)，其中參數(shù)估計(jì)是洪水頻率分析的重要部分，但未知參數(shù)的估計(jì)具有不確定性，結(jié)果的可靠性受到不確定性的影響，因此在參數(shù)估計(jì)中需要盡可能綜合各種現(xiàn)有信息來(lái)進(jìn)行統(tǒng)計(jì)推斷。貝葉斯統(tǒng)計(jì)學(xué)與經(jīng)典統(tǒng)計(jì)學(xué)的主要區(qū)別就是將先驗(yàn)信息加以利用，與樣本信息結(jié)合來(lái)進(jìn)行統(tǒng)計(jì)推斷，可以更為客觀全面地描述未知參數(shù)的統(tǒng)計(jì)特性。后驗(yàn)分布的推求是貝葉斯推斷的關(guān)鍵步驟，但常常因其中的高維積分太過(guò)復(fù)雜而難以求解。馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，簡(jiǎn)稱MCMC）方法是解決復(fù)雜積分問(wèn)題的有效途徑。本文選用3種MCMC抽樣算法，結(jié)合實(shí)例對(duì)不同算法在P-III分布參數(shù)估計(jì)中的計(jì)算效率進(jìn)行了比較和評(píng)估。

1 貝葉斯理論

基于總體信息、樣本信息和先驗(yàn)信息進(jìn)行的統(tǒng)計(jì)推斷稱為貝葉斯統(tǒng)計(jì)學(xué)，其與經(jīng)典統(tǒng)計(jì)學(xué)的主要區(qū)別為是否利用先驗(yàn)信息。同時(shí)，貝葉斯學(xué)派重視已出現(xiàn)的樣本觀察值，不考慮尚未發(fā)生的樣本觀察值。先驗(yàn)信息是在抽樣之前有關(guān)統(tǒng)計(jì)問(wèn)題的一些信息，主要來(lái)源于經(jīng)驗(yàn)和歷史資料，經(jīng)過(guò)挖掘、加工和數(shù)量化后形成的先驗(yàn)分布可用于后驗(yàn)分布的推求，基于后驗(yàn)分布對(duì)參數(shù)進(jìn)行統(tǒng)計(jì)推斷更為合理有效，因?yàn)楹篁?yàn)分布可以看作人們用總體信息和樣本信息對(duì)先驗(yàn)分布進(jìn)行調(diào)整的結(jié)果，反映了人們?cè)诔闃雍髮?duì)參數(shù)的認(rèn)識(shí)[1]。其基本觀點(diǎn)是：將未知參數(shù)視為隨機(jī)變量而非復(fù)雜的未知常量[2]。某站點(diǎn)洪水的密度函數(shù)f（x，θ）在貝葉斯統(tǒng)計(jì)中可記為f（x|θ），其中θ為未知參數(shù)矢量，表示在未知參數(shù)矢量給定某個(gè)值時(shí)，總體指標(biāo)X的條件分布，設(shè)先驗(yàn)分布為 π（θ），則后驗(yàn)分布 π（θ｜x）的表達(dá)式為：

2 MCMC方法

馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，簡(jiǎn)稱MCMC）方法[3]目前是處理復(fù)雜統(tǒng)計(jì)問(wèn)題的有效手段，尤其在經(jīng)常需要計(jì)算復(fù)雜高維積分的貝葉斯統(tǒng)計(jì)領(lǐng)域更是如此。在貝葉斯統(tǒng)計(jì)學(xué)中，高維積分主要用于求得普通方法無(wú)法獲得的后驗(yàn)分布密度。MCMC方法的基本思想是通過(guò)建立一個(gè)平穩(wěn)分布為π（θ）的馬爾科夫鏈來(lái)得到π（θ）的樣本，基于樣本來(lái)進(jìn)行各種統(tǒng)計(jì)推斷。若定義和實(shí)施合理，馬爾科夫鏈將會(huì)收斂，此時(shí)其極限分布即為后驗(yàn)分布。其步驟如下：

1）在觀測(cè)樣本上選擇一個(gè)平穩(wěn)分布為π（θ）的馬爾科夫鏈，使其轉(zhuǎn)移核為 f（*，*）；

2）由觀測(cè)樣本上的某一點(diǎn)θ0出發(fā)，用（1）中的馬爾科夫鏈產(chǎn)生點(diǎn)序列θ1，…，θn；

3）對(duì)某個(gè)m和足夠大的n，任一函數(shù)f（x）的期望估計(jì)如下：

可以看出采用MCMC方法時(shí)，構(gòu)造轉(zhuǎn)移核是至關(guān)重要的，不同的MCMC方法其實(shí)質(zhì)是轉(zhuǎn)移核的構(gòu)造方法不同。

2.1 MCMC抽樣算法

2.1.1 MH算法

Metropolis-Hastings（MH）抽樣算法是通過(guò)服從建議分布的候選點(diǎn)γ獲得t+1時(shí)刻的狀態(tài)。一般要求建議分布的形式盡可能簡(jiǎn)單，易于計(jì)算，馬爾科夫鏈下一個(gè)狀態(tài)被接受為候選點(diǎn)的概率為：

抽樣步驟如下：

1）給定初值 θ0，令 t=0；

3）從均勻分布［0，1］產(chǎn)生 U；

4）若 U≤α（Xt，γ），令 Xt+1=γ，否則 Xt+1=Xt；

5）令 t=t+1，重復(fù)步驟（2）～（5），直到產(chǎn)生足夠的抽樣樣本為止。

2.1.2 DR算法

延遲拒絕（Delayed Rejection，DR）抽樣算法中，目標(biāo)分布為π時(shí)，對(duì)密度為的候選值θ*，其迭代接受概率為：

當(dāng)θ*被拒絕時(shí)，下一個(gè)密度為q2（θ，θ*，θ**）的候選值θ**的迭代接受概率為：

2.1.3 AM算法

與MH算法和DR算法相比，自適應(yīng)Metropolis（Adaptive Metropolis，AM）抽樣算法并不需要指定參數(shù)的建議分布，而是由后驗(yàn)參數(shù)的協(xié)方差矩陣在每一次迭代后自適應(yīng)地調(diào)整來(lái)得到參數(shù)樣本[4,5]，其抽樣步驟如下：

1）初始化t=0，給定協(xié)方差矩陣C0；

2）在每一步迭代中產(chǎn)生新參數(shù)θ*；

4）產(chǎn)生[0，1]區(qū)間均勻分布隨機(jī)數(shù)r，確定接受新參數(shù)，即

5）重復(fù)步驟（2）～（4），直到產(chǎn)生足夠的抽樣樣本為止。

2.2 MCMC方法的收斂性診斷

MCMC方法的收斂性診斷方式有多種，方法一是對(duì)多個(gè)初值產(chǎn)生多個(gè)馬爾科夫鏈，觀察歷史迭代圖，波動(dòng)幅度穩(wěn)定時(shí)說(shuō)明收斂；方法二是觀察遍歷均值，在得到的馬爾科夫鏈中每隔一段距離計(jì)算參數(shù)的遍歷均值，當(dāng)這樣算得的均值穩(wěn)定時(shí)，可認(rèn)為抽樣收斂；方法三是方差比收斂性診斷，當(dāng)方差比曲線穩(wěn)定在1.0左右時(shí)說(shuō)明抽樣收斂。本文采用陸樂(lè)等[6]提出的EM（Ensemble Mean）和 ES（Ensemble Spread）兩個(gè)指標(biāo)來(lái)診斷樣本是否收斂，當(dāng)樣本序列的EM和ES值穩(wěn)定時(shí)認(rèn)為抽樣收斂。

式中，np為參數(shù)個(gè)數(shù)，Meani和Vari分別為第i個(gè)參數(shù)樣本序列的均值和方差。

3 P-III分布參數(shù)的貝葉斯估計(jì)

P-III分布的概率密度為：

式中，Γ（·）為 gamma函數(shù)，α、β、ɑ0分別為形狀、尺度和位置參數(shù)。它們與樣本均值x、變差系數(shù)Cv和偏態(tài)系數(shù)Cs的關(guān)系為：

后驗(yàn)分布為：

使用MH算法和DR算法抽樣時(shí)，選取多元正態(tài)分布為建議分布。

4 應(yīng)用實(shí)例

為比較不同MCMC抽樣算法的計(jì)算效率，選擇陜北地區(qū)劉家河、神木等5個(gè)水文測(cè)站的年最大洪峰流量系列資料進(jìn)行應(yīng)用，資料概況如表1所示。

設(shè)置模擬抽樣次數(shù)為15 000次，通過(guò)Matlab編制程序，采用MH、DR和AM三種抽樣算法分別進(jìn)行模擬抽樣，獲得模擬參數(shù)樣本系列。所得樣本序列的EM、ES曲線及經(jīng)驗(yàn)頻率和理論頻率擬合結(jié)果如圖1～圖5所示。

表1 陜北地區(qū)5個(gè)測(cè)站洪峰流量系列長(zhǎng)度

圖1 劉家河站EM-ES-累積概率比較

圖2 神木站EM-ES-累積概率比較

圖3 綏德站EM-ES-累積概率比較

圖4 趙石窯站EM-ES-累積概率比較

圖5 志丹站EM-ES-累積概率比較

由圖可知，當(dāng)產(chǎn)生15000個(gè)參數(shù)樣本后，劉家河站AM抽樣的ES值大約穩(wěn)定在2.5左右，且抽樣次數(shù)不足5 000次時(shí)曲線便已穩(wěn)定，DR抽樣的ES值在2.0附近波動(dòng)，而MH抽樣的ES曲線仍呈上升趨勢(shì)，說(shuō)明馬爾科夫鏈尚未收斂到參數(shù)的后驗(yàn)分布；EM曲線趨勢(shì)與ES曲線相同，其余各站情況類似，即在同等抽樣次數(shù)下，AM抽樣的ES值最早達(dá)到穩(wěn)定，收斂最快，DR抽樣的結(jié)果次之，MH抽樣在試驗(yàn)次數(shù)內(nèi)無(wú)法收斂到參數(shù)的后驗(yàn)分布。對(duì)于理論頻率和經(jīng)驗(yàn)頻率的擬合結(jié)果，由圖可知，三種抽樣方法的結(jié)果基本相同，可認(rèn)為不同抽樣算法對(duì)參數(shù)樣本的頻率擬合效果沒(méi)有影響。

5 結(jié)論

本文基于貝葉斯方法，比較了不同MCMC抽樣算法對(duì)P-III分布參數(shù)估計(jì)的效率，結(jié)論如下：

1）分別采用MH、DR、AM抽樣算法時(shí)，AM算法最快收斂，DR算法次之，MH算法收斂最慢，其中AM算法不依賴于參數(shù)的建議分布，能夠獲得較快收斂速度，且收斂指標(biāo)曲線最為平穩(wěn)，顯著優(yōu)于其余2種算法。

2）不同算法抽樣對(duì)所得參數(shù)樣本經(jīng)驗(yàn)頻率和理論頻率的擬合效果沒(méi)有影響。