歐陽(yáng)資生,甘 柳
在各種水利規(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í)證分析。
通過(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ì)α?。
現(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)意的。
作為模型的一個(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í)間序列圖
根據(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ì)圖
以上四個(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.