洪水頻率分布的尾指數估計

歐陽資生，甘 柳

0 引言

在各種水利規劃設計中，都需要進行水文頻率分析工作，以保證達到工程要求的設計值。而在洪水頻率分析中，無非是討論像洪水或超出某一警戒水平的水位或流量的分布情況。在極值統計中，像洪水等這類事件發生頻率很低，但是一經發生就容易造成較大損失的事件叫極值事件。事實上，在洪災風險管理中，如何發現這些極值事件的發生概率和某個極值分位數點對風險管理者是相當重要的。我們知道，從極值理論的角度來看，這些極值事件的概率和極值分位數被一種稱為極值分布的尾指數所控制。對于正態分布，它的尾部呈指數函數衰減，其尾指數為零。當尾指數大于零時，分布尾部呈冪函數衰減，我們稱其為厚尾分布。厚尾模型在諸如金融、保險、水利等很多場合都是一個應用非常廣泛的分布模型。厚尾分布中，尾指數越大，其尾部越厚。因此尾指數可作為衡量洪水分布規律的重要指標。

假設X1,X2,???,Xn是一列正的，獨立同分布的隨機變量序列，具有共同的分布函數F(x)：

則分布函數F(x)被稱為厚尾的。模型（1）的F(x)也稱為Pareto型分布。?F是無窮遠處的緩變函數，滿足對所有的η＞0：

這里，γ稱為極值指數，而如果令α=1/γ，則α稱為尾指數。尾指數或極值指數的估計目前仍然是極值統計的一個熱點問題。如Danielsson(2001)利用自助法研究了尾指數的估計問題。Beirlant et al.(2008)提出了一種修正的偏差減少方法對尾指數和極值分位數進行估計。Beran and Schell(2010)在小樣本情況下構造了一個穩健的M估計方法對尾指數進行了估計。Brito and Freitas(2010)研究了相依數據尾指數估計的相合性問題。但總的說來，對尾指數或極值指數的的估計，不外乎是Pickands估計、Hill估計和矩估計三種估計及其拓展估計方法。在這三種估計中，Hill估計和矩估計實際應用中相對較多，而Pickands估計實用性并不強，基本不用。

如果設：

為 n個觀測值 X1,X2,???,Xn的順序統計量，Hill估計定義為：

Hill估計雖然在理論上具有很好的大樣本性質，但是在實際應用中，并不好操作，我們可以從圖1就可很容易理解。

圖1 極值指數的Hill估計圖

圖1 是我們利用學生-t4分布對Hill估計結果作的一個隨機模擬，我們作了200次隨機模擬，每次模擬的樣本量是500，我們給出了模擬的200次估計的1/4分位數，中位數和3/4分位數，在學生-t4分布中，極值指數的真值γ=0.25，但是極值指數的真值到底取多少，我們從圖1中很難做出判斷。換句話說，我們不知道門限值取多大時才能對樣本進行有效分割。事實上，如何選取合適的門限值是估計極值指數或尾指數的基礎，也是我們進行洪水頻率分析的必要程序。

本文中，我們將基于指數回歸模型，給出矩估計的門限值和樣本點分割的選取原理和方法，然后利用MC方法進行模擬說明門限值選取的合理性，最后利用所構建的模型對湖南省四個水文觀測站的水文數據進行實證分析。

1 尾指數估計的統計建模

通過前面的模擬，我們看到Hill估計并不好操作，門限值選擇不當將導致極值指數估計的較大偏差。事實上，極值指數三種估計方法中，矩估計相對來說較為穩健，實際應用中也相對較多。因此，本文我們采用矩估計來估計極值指數和尾指數。

矩估計由Dekkers、Einmahl&de Hall(1989)提出，其定義為：對 k∈{3,???,(n-1)}

盡管矩估計較Hill估計穩健，在理論上也具有更好的大樣本性質，但是在實際應用中，也同樣存在門限值選取的問題。我們可以從圖2就很容易發現，雖然相對Hill估計而言，矩估計要穩健得多，但是門限值的選取問題仍然是一個需要解決的問題。

圖2 極值指數的矩估計圖

圖2 是我們利用學生-t1分布對矩估計結果作的一個隨機模擬，我們同樣作了200次隨機模擬，每次模擬的樣本量大小為500，我們給出了模擬的500次估計的中位數。在學生-t1分布中，極值指數的真值為γ=1。雖然比起Hill估計而言要穩健得多，但是，極值指數的真值到底取多少，我們從圖2中也不好判斷。換句話說，我們不知道門限值取多少時，才能對樣本實現最優分割。

和其他估計一樣，在用矩估計對尾指數進行估計時，首先是確定門限值，找出超出門限值以上的觀察數據；也就是對所觀察到的樣本值的順序統計量進行有效分割，得到用于估計的觀察數據，然后才能進行估計。但需要指出的是，門限值的選取問題卻一直是困擾極值工作者的一個難題。門限值越大，可以分析的數據越少，這時，被分析的數據比較接近分布的極端，分析的偏差減少，但由于數據過少，估計的方差增加；反之，門限值過小，被分析的數據增加，分析的方差減少，但偏差卻增加了。對這個問題的研究，統計工作者提出了許多方案。如Dupuis(1998)建議從參數的穩健性出發來確定門限值；Guillou(2001)、Matthys&Beirlant(2003),Beirlant et al.(1996,2004)，歐陽資生(2008)等建議使用最小化均方誤差或漸近二階矩來獲得門限值；Gomes et al.(2008)建議使二階參數估計的偏差達到最小從而通過一個啟發式適應過程得到門限選擇方法；Vandewalle et al.(2008)通過使PDC估計(partial density component estimation)的積分均方誤差達到最小來獲得門限值，對樣本進行分割。

采用SPSS 19.0軟件對數據進行分析處理，計量資料以（均數±標準差）表示，采用t檢驗；計數資料以（n，%）表示，采用χ2檢驗，以P＜0.05表示差異具有統計學意義。

下面：我們將基于指數回歸模型，在漸近最小均方誤差的準則下，給出矩估計的門限值和樣本點分割的選取原理和方法，并提出極值指數和尾指數估計的算法。

對于隨機變量序列X1,X2,???,Xn的分布函數F(x)，如前所述，我們假設F(x)是Pareto型的，其原因主要是基于Pareto型在極值分布中的地位和作用。我們知道，這種分布在金融、保險、水利中都被廣泛應用。例如，在巨災統計數據中，有一個廣為人知的事實，即巨災統計數據是厚尾的，因此，可以直接假設巨災統計數據分布服從Pareto分布。

對于模型(1)中的緩變函數，有一個被廣泛接受的假設：

假設?λ:存在一個實常數ρ＜0和一個正的比率函數b(x)，滿足當x→∞時，b(x)→0，且使得對所有的η≥1,

其中 kρ(η)= ∫1ηvρ-1d v=(ηρ-1)/ρ ，若 ρ=0 ，則 kρ(η)=log(η)。需要說明的是，假設?λ條件并不苛刻，一般的緩變函數均能滿足這個條件。

在假設?λ下，我們按照Beirlant et al.(2004)，歐陽資生(2008)的指數回歸模型方法來選取k，從而進一步確定參數 γ?k,b?n,k,ρ?k。為此，建立如下指數回歸模型：

這里，f1,f2,???,fk是一列獨立的，服從標準指數分布的隨機變量。在式(6)中，利用最大似然估計，得到參數γ,bn,k,ρ的估計值：

類似于Beirlant et al.(2004)，歐陽資生(2008)，我們可得在指數回歸模型中，極值指數用矩估計作為估計量時的AMSE為：

因此，樣本的最優分割k?optn為：

因此，根據以上原理，我們可得基于指數回歸模型的樣本分割方法，進而得到尾指數的矩估計的算法，算法如下：

(1)對指數回歸模型式(6)，利用極大似然估計，對k∈{3,???,(n-1)}計算參數 γ ,bn,k,ρ 的估計值{(γ?k,b?n,k,ρ?k),k∈{3,???,(n-1)}

(2)對 k∈{3,???,(n-1)}計算 AMSE(γ?Mk)

(3)利用

獲得 k?optn

(4)根據矩估計式(4)和步驟3的最優k?optn，可得極值指數的最優估計和尾指數估計α?。

2 隨機模擬

現在，為驗證我們的模型，我們對下列極值分布進行蒙特卡洛模擬：

(1)Burr(1,1,1)分布。 Burr(θ,τ,λ)分布的分布函數滿足：

(2)Burr(1,0.5,2)分布

(3)Frechet(1)分布.Fre ch et(γ)分布的分布函數滿足：

(4)Frechet(2)分布

(5)學生-t4分布

表1 極值指數估計及其誤差估計模擬結果表

在蒙特卡洛模擬時，我們對每一種分布作了500次模擬，每次模擬的樣本量均為1000。表1分別給出了500次模擬中相應的最優k值、γ的估計值的平均及其標準差、AMSE的平均。從表1可以看出，在矩估計中，借助于指數回歸模型獲得門限值、樣本點分割方法和極值指數估計值，其結果是令人非常滿意的。

3 湖南省洪水頻率分析的實證分析

3.1 數據描述

作為模型的一個應用，我們對洞庭湖周邊的桃源、津市、沙頭、石龜山等四個站點的水文數據中的水流量進行實證分析。數據跨度為1998年元月1日至2010年4月1日共4316個日數據。為對數據的基本情況有一直觀了解，我們在表2中列出了相應的統計量。同時，也繪畫了其相應的時間序列圖(圖3)。從表2可以發現，這四個站點的數據均呈現明顯的厚尾現象，同時，從圖3也可看出其波動明顯。

表2 洞庭湖周邊的四個水文觀測站水流量基本統計特征

圖3 （1998.1.1-2010.4.1）洞庭湖周邊四個水文觀測站水流量時間序列圖

3.2 尾指數的估計

根據前文的極值指數估計模型和計算方法，我們首先可以得到樣本的最優分割方法從而得到超出門限值的樣本個數，然后得到極值指數的估計值，最后利用α?=1/γ?即可得湖南省四個水文觀測站點的水流量分布的尾指數估計。

在圖4中，我們給出了k=1,2,...,4310時極值指數的估計圖,其中虛線代表了最優的k值。而表3分別給出了這四個站點在最優的k值下的極值指數和尾指數最優估計值，從表3可以看出，桃源站的尾指數是最大的，津市的尾指數是最小的。

兩點說明：(1)從前文對學生-t1分布的極值指數的隨機模擬圖2可看出，即使是估計時相對穩健的矩估計也不好直接判斷尾指數的真值，這點在圖4中也得到了印證。從圖4中，如果沒有一定的準則，我們是無法獲知在何時對樣本實施分割，也就無法得到各個站點水流量數據的尾指數的真值了。但是，如果我們借助于指數回歸模型，采用使矩的AMSE達到最小作為評價標準，我們就可以很好地解決一個問題。(2)正如Beirlant et al.(1996,2004)在借助指數回歸模型對尾指數進行Hill估計時作的評述，指數回歸模型相對來說，較為穩健。同時，由于有效的利用了極大似然估計，因此計算的速度也較快，這也是我們在進行矩估計時，借助指數回歸模型進行建模的主要原因。

圖4 洞庭湖周邊的桃源等四個水文站點水流量數據的極值指數估計圖

4 結論

以上四個水文觀測站點都是位于洞庭湖地區周邊站點，由前面的分析結果可以發現，四水系流域的水位變化情況均為厚尾分布，都可以通過極值分布加以較好地擬合。當然，流經不同站點的水流量是不一樣的，且不同規模洪水流量的變化幅度亦有所區別，因此在實施防洪措施時應實事求是，依據不同的情況有區別的對待，這樣才能既做到全面有效防洪減災又能盡可能的降低不必要投入，減少浪費。

[1]Beirlant,J.,Figueiredo,F.,Gomes,M.I.,Vandewalle,B.Improved Re?duced-bias Tail Index and Quantile Estimators[J].J.Statist.Plann.and Inference，2008,138.

[2]Beirlant,J.,Goegebeur,Y.,Segers,J.,Teugels,J.Statistics of Ex?tremes.Theory and Applications[M].NewYork:Wiley,2004.

[3]Beirlant,J.,Vynckier,P.,Teugels,J.L.Tail Index Estimation,Pareto Quantile Plots,and Regression Diagnostics[J].J.Amer.Statist.Assoc，1996,91.

[4]Beran,J.,Schell,D.On Robust Tail Index Estimation.Computational Statisticsand Data Analysis,doi:10.1016/j.csda[J].2010.

[5]Brito,M.,Freitas,A.C.Consistent Erstimation of the Tail Index for De?pendent Data[J].Statistics and Probability Letters，2010,(80).

[6]Danielsson,J.Using a Bootstrap Method to Choose the Sample Fraction in Tail Index Estimation[J].Journal of Multivariate Analysis，2001,76.

[7]Dekkers,A.,de Haans,L.A Moment Estimator for the Index of an Ex?treme-value Distribution[J].Ann Statist,1989,17(4).

[8]Dupuis,D.J.Exceedances over High Thresholds:a Guide to Thresh?old Selection[J].Extremes,1998,3(1).

[9]Gomes,M.I.,Henriques Rodrigues,L,Vandewalle,B.，Viseu,C.A Heu?ristic Adaptive Choice of the Threshold for Bias-corrected Hill Esti?mators[J].J.Statist.Comput.And Simulation,2008,78(2).

[10]Guillou,A.,Hall,P.A Diagnostic for Selecting the Threshold in Ex?treme Analysis[J].J.R.Statist.Soc.Ser B,2001,63.

[11]Matthys,G.,Beirlant,J.Estimating the Extreme Value Index and High Quantiles with Exponential Regression Models[J].Statistica Si?nica,2003,13.

[12]Vandewalle,B.,Beirlant,J.,Christmann,A.,Hubert,M.A Robust Estimator for the Tail Index of Pareto-type Distributions[J].Compu?tational Statistics&Data Analysis,2007,51.

[13]歐陽資生.厚尾分布的極值分位數估計與極值風險測度研究[J].數理統計與管理,2008,27.