歐陽(yáng)資生，甘 柳

0 引言

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

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

則分布函數(shù)F(x)被稱(chēng)為厚尾的。模型（1）的F(x)也稱(chēng)為Pareto型分布。?F是無(wú)窮遠(yuǎn)處的緩變函數(shù)，滿(mǎn)足對(duì)所有的η＞0：

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

如果設(shè)：

為 n個(gè)觀測(cè)值 X1,X2,???,Xn的順序統(tǒng)計(jì)量，Hill估計(jì)定義為：

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

圖1 極值指數(shù)的Hill估計(jì)圖

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

本文中，我們將基于指數(shù)回歸模型，給出矩估計(jì)的門(mén)限值和樣本點(diǎn)分割的選取原理和方法，然后利用MC方法進(jìn)行模擬說(shuō)明門(mén)限值選取的合理性，最后利用所構(gòu)建的模型對(duì)湖南省四個(gè)水文觀測(cè)站的水文數(shù)據(jù)進(jìn)行實(shí)證分析。

1 尾指數(shù)估計(jì)的統(tǒng)計(jì)建模

通過(guò)前面的模擬，我們看到Hill估計(jì)并不好操作，門(mén)限值選擇不當(dāng)將導(dǎo)致極值指數(shù)估計(jì)的較大偏差。事實(shí)上，極值指數(shù)三種估計(jì)方法中，矩估計(jì)相對(duì)來(lái)說(shuō)較為穩(wěn)健，實(shí)際應(yīng)用中也相對(duì)較多。因此，本文我們采用矩估計(jì)來(lái)估計(jì)極值指數(shù)和尾指數(shù)。

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

盡管矩估計(jì)較Hill估計(jì)穩(wěn)健，在理論上也具有更好的大樣本性質(zhì)，但是在實(shí)際應(yīng)用中，也同樣存在門(mén)限值選取的問(wèn)題。我們可以從圖2就很容易發(fā)現(xiàn)，雖然相對(duì)Hill估計(jì)而言，矩估計(jì)要穩(wěn)健得多，但是門(mén)限值的選取問(wèn)題仍然是一個(gè)需要解決的問(wèn)題。

圖2 極值指數(shù)的矩估計(jì)圖

圖2 是我們利用學(xué)生-t1分布對(duì)矩估計(jì)結(jié)果作的一個(gè)隨機(jī)模擬，我們同樣作了200次隨機(jī)模擬，每次模擬的樣本量大小為500，我們給出了模擬的500次估計(jì)的中位數(shù)。在學(xué)生-t1分布中，極值指數(shù)的真值為γ=1。雖然比起Hill估計(jì)而言要穩(wěn)健得多，但是，極值指數(shù)的真值到底取多少，我們從圖2中也不好判斷。換句話(huà)說(shuō)，我們不知道門(mén)限值取多少時(shí)，才能對(duì)樣本實(shí)現(xiàn)最優(yōu)分割。

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

采用SPSS 19.0軟件對(duì)數(shù)據(jù)進(jìn)行分析處理，計(jì)量資料以（均數(shù)±標(biāo)準(zhǔn)差）表示，采用t檢驗(yàn)；計(jì)數(shù)資料以（n，%）表示，采用χ2檢驗(yàn)，以P＜0.05表示差異具有統(tǒng)計(jì)學(xué)意義。

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

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

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

假設(shè)?λ:存在一個(gè)實(shí)常數(shù)ρ＜0和一個(gè)正的比率函數(shù)b(x)，滿(mǎn)足當(dāng)x→∞時(shí)，b(x)→0，且使得對(duì)所有的η≥1,

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

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

這里，f1,f2,???,fk是一列獨(dú)立的，服從標(biāo)準(zhǔn)指數(shù)分布的隨機(jī)變量。在式(6)中，利用最大似然估計(jì)，得到參數(shù)γ,bn,k,ρ的估計(jì)值：

類(lèi)似于Beirlant et al.(2004)，歐陽(yáng)資生(2008)，我們可得在指數(shù)回歸模型中，極值指數(shù)用矩估計(jì)作為估計(jì)量時(shí)的AMSE為：

因此，樣本的最優(yōu)分割k?optn為：

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

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

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

(3)利用

獲得 k?optn

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

2 隨機(jī)模擬

現(xiàn)在，為驗(yàn)證我們的模型，我們對(duì)下列極值分布進(jìn)行蒙特卡洛模擬：

(1)Burr(1,1,1)分布。 Burr(θ,τ,λ)分布的分布函數(shù)滿(mǎn)足：

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

(3)Frechet(1)分布.Fre ch et(γ)分布的分布函數(shù)滿(mǎn)足：

(4)Frechet(2)分布

(5)學(xué)生-t4分布

表1 極值指數(shù)估計(jì)及其誤差估計(jì)模擬結(jié)果表

在蒙特卡洛模擬時(shí)，我們對(duì)每一種分布作了500次模擬，每次模擬的樣本量均為1000。表1分別給出了500次模擬中相應(yīng)的最優(yōu)k值、γ的估計(jì)值的平均及其標(biāo)準(zhǔn)差、AMSE的平均。從表1可以看出，在矩估計(jì)中，借助于指數(shù)回歸模型獲得門(mén)限值、樣本點(diǎn)分割方法和極值指數(shù)估計(jì)值，其結(jié)果是令人非常滿(mǎn)意的。

3 湖南省洪水頻率分析的實(shí)證分析

3.1 數(shù)據(jù)描述

作為模型的一個(gè)應(yīng)用，我們對(duì)洞庭湖周邊的桃源、津市、沙頭、石龜山等四個(gè)站點(diǎn)的水文數(shù)據(jù)中的水流量進(jìn)行實(shí)證分析。數(shù)據(jù)跨度為1998年元月1日至2010年4月1日共4316個(gè)日數(shù)據(jù)。為對(duì)數(shù)據(jù)的基本情況有一直觀了解，我們?cè)诒?中列出了相應(yīng)的統(tǒng)計(jì)量。同時(shí)，也繪畫(huà)了其相應(yīng)的時(shí)間序列圖(圖3)。從表2可以發(fā)現(xiàn)，這四個(gè)站點(diǎn)的數(shù)據(jù)均呈現(xiàn)明顯的厚尾現(xiàn)象，同時(shí)，從圖3也可看出其波動(dòng)明顯。

表2 洞庭湖周邊的四個(gè)水文觀測(cè)站水流量基本統(tǒng)計(jì)特征

圖3 （1998.1.1-2010.4.1）洞庭湖周邊四個(gè)水文觀測(cè)站水流量時(shí)間序列圖

3.2 尾指數(shù)的估計(jì)

根據(jù)前文的極值指數(shù)估計(jì)模型和計(jì)算方法，我們首先可以得到樣本的最優(yōu)分割方法從而得到超出門(mén)限值的樣本個(gè)數(shù)，然后得到極值指數(shù)的估計(jì)值，最后利用α?=1/γ?即可得湖南省四個(gè)水文觀測(cè)站點(diǎn)的水流量分布的尾指數(shù)估計(jì)。

在圖4中，我們給出了k=1,2,...,4310時(shí)極值指數(shù)的估計(jì)圖,其中虛線(xiàn)代表了最優(yōu)的k值。而表3分別給出了這四個(gè)站點(diǎn)在最優(yōu)的k值下的極值指數(shù)和尾指數(shù)最優(yōu)估計(jì)值，從表3可以看出，桃源站的尾指數(shù)是最大的，津市的尾指數(shù)是最小的。

兩點(diǎn)說(shuō)明：(1)從前文對(duì)學(xué)生-t1分布的極值指數(shù)的隨機(jī)模擬圖2可看出，即使是估計(jì)時(shí)相對(duì)穩(wěn)健的矩估計(jì)也不好直接判斷尾指數(shù)的真值，這點(diǎn)在圖4中也得到了印證。從圖4中，如果沒(méi)有一定的準(zhǔn)則，我們是無(wú)法獲知在何時(shí)對(duì)樣本實(shí)施分割，也就無(wú)法得到各個(gè)站點(diǎn)水流量數(shù)據(jù)的尾指數(shù)的真值了。但是，如果我們借助于指數(shù)回歸模型，采用使矩的AMSE達(dá)到最小作為評(píng)價(jià)標(biāo)準(zhǔn)，我們就可以很好地解決一個(gè)問(wèn)題。(2)正如Beirlant et al.(1996,2004)在借助指數(shù)回歸模型對(duì)尾指數(shù)進(jìn)行Hill估計(jì)時(shí)作的評(píng)述，指數(shù)回歸模型相對(duì)來(lái)說(shuō)，較為穩(wěn)健。同時(shí)，由于有效的利用了極大似然估計(jì)，因此計(jì)算的速度也較快，這也是我們?cè)谶M(jìn)行矩估計(jì)時(shí)，借助指數(shù)回歸模型進(jìn)行建模的主要原因。

圖4 洞庭湖周邊的桃源等四個(gè)水文站點(diǎn)水流量數(shù)據(jù)的極值指數(shù)估計(jì)圖

4 結(jié)論

以上四個(gè)水文觀測(cè)站點(diǎn)都是位于洞庭湖地區(qū)周邊站點(diǎn)，由前面的分析結(jié)果可以發(fā)現(xiàn)，四水系流域的水位變化情況均為厚尾分布，都可以通過(guò)極值分布加以較好地?cái)M合。當(dāng)然，流經(jīng)不同站點(diǎn)的水流量是不一樣的，且不同規(guī)模洪水流量的變化幅度亦有所區(qū)別，因此在實(shí)施防洪措施時(shí)應(yīng)實(shí)事求是，依據(jù)不同的情況有區(qū)別的對(duì)待，這樣才能既做到全面有效防洪減災(zāi)又能盡可能的降低不必要投入，減少浪費(fèi)。

[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]歐陽(yáng)資生.厚尾分布的極值分位數(shù)估計(jì)與極值風(fēng)險(xiǎn)測(cè)度研究[J].數(shù)理統(tǒng)計(jì)與管理,2008,27.