程月,田茂再,2

（1.新疆財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與數(shù)據(jù)科學(xué)學(xué)院,新疆 烏魯木齊 830012；2.中國人民大學(xué)應(yīng)用統(tǒng)計(jì)科學(xué)研究中心,北京 100190）

統(tǒng)計(jì)分布被廣泛用于分析和預(yù)測真實(shí)世界的現(xiàn)象,但由于新的應(yīng)用程序和事件的快速出現(xiàn),需要不斷開發(fā)新的分布,因此統(tǒng)計(jì)分布理論和產(chǎn)生新的分布引起了人們極大的研究興趣.近幾年,研究者對極端事件的極端值在各個(gè)領(lǐng)域進(jìn)行分析,尤其是在數(shù)學(xué)、醫(yī)學(xué)以及統(tǒng)計(jì)學(xué)等領(lǐng)域.雖然Gumbel Type-II 分布已成為廣義極值理論的重要組成部分,但與Weibull 分布不同的是該分布由于擬合欠佳而導(dǎo)致其應(yīng)用的普及率較低.由于在現(xiàn)實(shí)生活實(shí)踐中經(jīng)常遇到較多復(fù)雜現(xiàn)象的危險(xiǎn)率是非常單調(diào)的,因此不能用Gumbel Type-II 分布進(jìn)行建模.相反,此分布僅應(yīng)用于具有單調(diào)失效率的實(shí)際數(shù)據(jù).

為了解決這一限制,目前國內(nèi)外有很多關(guān)于Gumbel Type-II（GT-II）分布最新發(fā)展和應(yīng)用的文獻(xiàn).Okorie 等[1]發(fā)現(xiàn)了具有三個(gè)參數(shù)的EG Type-2 分布,并利用一組實(shí)際數(shù)據(jù)展示了EG Type-2 分布在相比較的分布中具有最佳擬合；Ogunde 等[2]提出了一種新的四個(gè)參數(shù)的EGTT 分布并推出密度、次序統(tǒng)計(jì)量和矩目函數(shù)等表達(dá)式,最后采用豬的癌癥緩解時(shí)間和生存時(shí)間兩組生命數(shù)據(jù)證實(shí)EGTT 分布在建模真實(shí)生命數(shù)據(jù)中的適用性,進(jìn)而發(fā)現(xiàn)EGTT 分布更加靈活；Sindhu 等[3]推廣了一種GLET-GTII 分布以及一些特定的統(tǒng)計(jì)屬性, 提出了運(yùn)用貝葉斯估計(jì)和極大似然估計(jì)來討論GLET-GTII 分布的參數(shù)估計(jì),MCMC 模擬分析結(jié)果與估計(jì)量是一致的.最后收集了中國1 月23 日-3 月28 日和歐洲3 月1 日-3 月30 日因COVID-19 的死亡人數(shù)數(shù)據(jù)分析GLET-GTII 分布的適用性,使用經(jīng)典的擬合優(yōu)度指標(biāo)評估GLET-GTII 分布優(yōu)于絕大數(shù)已知的分布；袁丹華等[4]討論了在逐次II 型截尾樣本下的Gumbel 極值分布,分別利用Bayes 估計(jì)和IFM法對該分布進(jìn)行參數(shù)估計(jì),結(jié)果表明當(dāng)參數(shù)偏離1 的程度較小時(shí),這兩種方法的估計(jì)出來的結(jié)果相差不大,但當(dāng)參數(shù)偏離程度較大時(shí),IFM方法的估計(jì)結(jié)果較差于Bayes 估計(jì)的結(jié)果；Sindhu 等[5]提出了一種具有三參數(shù)的ETGT-II 分布并分析了具體的統(tǒng)計(jì)屬性以及相關(guān)測度、矩生成函數(shù)、生存函數(shù)、分位數(shù)函數(shù)等.最后為了檢驗(yàn)ETGT-II 分布的有效性進(jìn)行仿真分析,利用已有文獻(xiàn)中收集的歐洲和中國因新冠肺炎死亡的人數(shù)數(shù)據(jù)驗(yàn)證了ETGT-II 分布相比其他分布具有較好的適用性；Willayat 等[6]提出了一種新的概率MOEGT-II 分布,它可以對各種形狀的失效率函數(shù)進(jìn)行建模,并利用MLE 方法對MOEGT-II 分布的未知參數(shù)進(jìn)行估計(jì),運(yùn)用3 種不同的迭代方法對不同大小樣本進(jìn)行模擬研究,發(fā)現(xiàn)與估計(jì)量是一致的.最后選取三組實(shí)際數(shù)據(jù)集與GT-II、EGT-II、MOGT-II、Frechet 等分布進(jìn)行比較分析,發(fā)現(xiàn)MOEGT-II 分布優(yōu)于其他分布且更具有適用性；Lone 等[7]使用了一種新的概率分布族與Gumbel Type-II 分布建立了一個(gè)獨(dú)特的NAPGT-II 分布,運(yùn)用兩組真實(shí)生活數(shù)據(jù)證明MAPGT-II 分布的P值最大,K-S 最低,比所對比的競爭分布產(chǎn)生了更好的擬合.本文主要使用MOAP 分布族與Gumbel Type-II 分布建立了一個(gè)新的分布,把這個(gè)新的分布命名為MOAPGT-II 分布.

1 理論方法

Marshall-Olkin 的分布家族是由Marshall 等[8]提出的一種新的方法,此方法增加一個(gè)形狀參數(shù)到分布中.假設(shè)G（x）和g（x）作為給定隨機(jī)變量X的累積分布函數(shù)（CDF）和概率密度函數(shù)（PDF）.則Marshall-Olkin家族的累積分布函數(shù)和概率密度函數(shù)分別為

Mahdavi 等[9]通過增加一個(gè)新參數(shù)α,以累積分布函數(shù)作為基線進(jìn)行變換,從而獲得更加靈活的分布族,由此提出了alpha power transform（APT）類的方法.假設(shè)G（x）是任何分布的累積分布函數(shù),則APT類的累積分布函數(shù)和概率密度函數(shù)分別為

Nassar 等[10]在前面兩位學(xué)者研究結(jié)果的基礎(chǔ)上,提出了一種新的方法,將額外的參數(shù)引入到分布族中,以獲得更大的靈活性,進(jìn)而得到一個(gè)新的分布,定義為Marshall Olkin alpha power（MOAP）分布族,表示MOAP 分布族的累積分布函數(shù)和概率密度函數(shù)分別為

本文假定基線分布為參數(shù)γ 和δ 的Gumbel Type-II 分布的累積分布函數(shù)（CDF）和概率密度函數(shù)（PDF）如下

通過使用MOAP 分布族和Gumbel Type-II 分布生成具有4 參數(shù)的MOAPGT-II 分布,其中且≠1,則MOAPGT-II 分布的累積分布函數(shù)（CDF）和概率密度函數(shù)（PDF）如下

不同參數(shù)值下的MOAPGT-II 分布的概率密度圖見圖1.

圖1 不同參數(shù)值下的MOAPGT-II 分布的概率密度圖Fig.1 Probability density diagram of MOAPGT-II distribution with different parameter values

在圖1 中可以得到,在其他參數(shù)不變的情況下,隨著參數(shù)α 的增大,控制了MOAPGT-II 分布函數(shù)從陡峭變?yōu)橹饾u平緩的變化.

首先,生存函數(shù)又叫可靠性函數(shù),它可以表示為一系列事件的隨機(jī)變量函數(shù),一般來說是用于表示一些基于時(shí)間的系統(tǒng)失敗或死亡概率.采用公式來確定MOAPGT-II 分布的生存函數(shù)如下

不同參數(shù)值下的MOAPGT-II 分布的生存函數(shù)見圖2.

圖2 不同參數(shù)值下的MOAPGT-II 分布的生存函數(shù)圖Fig.2 The survival function diagram of MOAPGT-II distribution with different parameter values

其次,通過以下公式來確定MOAPGT-II 分布的風(fēng)險(xiǎn)函數(shù)為

不同參數(shù)值下的MOAPGT-II 分布的風(fēng)險(xiǎn)函數(shù)見圖3.

圖3 不同參數(shù)值下的MOAPGT-II 分布的風(fēng)險(xiǎn)函數(shù)圖Fig.3 The risk function diagram of MOAPGT-II distribution under different parameter values

緊接著通過公式（13）來確定MOAPGT-II 分布的反風(fēng)險(xiǎn)函數(shù)為

最后通過公式（14）來確定MOAPGT-II 分布的賠率比為

2 極大似然估計(jì)

在本節(jié)中,討論使用R 語言編程來實(shí)現(xiàn)極大似然估計(jì)方法對所提出的MOAPGT-II 分布的參數(shù)進(jìn)行估計(jì).設(shè)y1,y2…,yn是樣本總體Y1,Y2…,Yn的一個(gè)樣本值,那么它的概率密度函數(shù)可以表示為f(yi,Ψ ),其中 Ψ為待估參數(shù),可以計(jì)算出聯(lián)合概率密度函數(shù)f(yi,Ψ )的似然函數(shù)為

對數(shù)似然函數(shù)的一般形式表示為

對（15）關(guān)于 Ψ求導(dǎo)有

lnL(Ψ )與L(Ψ )可以在同一點(diǎn)處取得極大值,因此 Ψ的極大似然估計(jì)值也可以從式（18）解得

采用極大似然方法對該分布的參數(shù)進(jìn)行估計(jì)的主要步驟如下：

（1）由總體分布導(dǎo)出樣本的聯(lián)合概率密度函數(shù)

（2）把樣本的聯(lián)合概率密度函數(shù)中的自變量看成已知的參數(shù),把待估參數(shù) Ψ看作自變量,從而得到對數(shù)似然函數(shù)lnL(Ψ)

（3）求解出對數(shù)似然函數(shù)lnL(Ψ )的最大值點(diǎn)的 Ψ值

（4）將樣本數(shù)據(jù)代入到式（19）的結(jié)果中,可以計(jì)算出參數(shù)的估計(jì)量.

一般來說,求解上面的非線性似然方程很難得到解析式,為了在數(shù)值上解決這些問題,需要采用牛頓法等迭代過程.因此本文使用optim（）R 函數(shù)執(zhí)行,參數(shù)方法=“L-BFGSB”.執(zhí)行后,分別以數(shù)字形式和視覺形式展示所獲得的結(jié)果.

MOAPGT-II 分布的模擬是通過選擇隨機(jī)樣本來實(shí)現(xiàn)的,利用模擬研究了MLE 相對于樣本大小n的性能,遵循以下步驟進(jìn)行數(shù)值模擬分析：設(shè)Yi,i=1,2,… ,n是來自MOAPGT-II分布模型式（10）中的獨(dú)立分布的樣本,估計(jì)程序生成的樣本量n 分別為50、100、200、400、600、800和1000,這些樣本是通過逆分布函數(shù)獲得的,設(shè)定整個(gè)估計(jì)程序的模擬次數(shù)m=1 000 次.最后考慮使用平均偏差（Bias）和均方誤差（MSE）作為評價(jià)標(biāo)準(zhǔn),其中是Ψ的估計(jì)值,則二者公式分別表示為

表1 中展示了MOAPGT-II 分布參數(shù)α、δ、θ和γ的初始值分別設(shè)置為1.05、0.1、1.6、4 的情況下的估計(jì)結(jié)果.

表1 數(shù)據(jù)集在不同分布情況下的參數(shù)估計(jì)和標(biāo)準(zhǔn)誤差Tab.1 Parameter estimation and standard error of data set in different distributions

由表1 可以看出,隨著樣本量的逐漸增大,各參數(shù)的平均值越來越接近于真實(shí)值,并趨于穩(wěn)定；隨著樣本量的增大,平均偏差逐漸減少并趨近于0；隨著樣本量的增大,均方誤差也逐漸接近于0,說明極大似然估計(jì)具有良好的穩(wěn)定性,表明本文用極大似然估計(jì)所提出的MOAPGT-II 分布的參數(shù)有效且具有一定的實(shí)際應(yīng)用價(jià)值.

3 應(yīng)用

在本節(jié)中,應(yīng)用一個(gè)真實(shí)的數(shù)據(jù)集來說明所提出的MOAPGT-II 分布的有效性,利用R 軟件中的MAXlik 包進(jìn)行驗(yàn)證.出于比較的目的,考慮采用以下幾種分布：①指數(shù)變換的Gumbel Type-II 分布（ETGT-II）、②廣義逆Weibull 分布（GIWD）、③Frechet 分布（Frechet）、④II 型耿貝爾分布（GT-II）、⑤指數(shù)Burr XII（EB-XII）分布.同時(shí)還參考了Eghwerido 等[11]基于MLE 估計(jì)的對數(shù)似然函數(shù)的不同分離方法：赤池信息量準(zhǔn)則（AIC）、改正的赤池信息量準(zhǔn)則（CAIC）和Hannan-Quin 信息量準(zhǔn)則（HQIC）.其中n 為樣本數(shù)量,k 代表待估計(jì)參數(shù)的數(shù)量,是取對數(shù)后的似然函數(shù)值.通常定義最優(yōu)的分布是給出以上準(zhǔn)則最小值的分布

本文選取了Sindhu 等[3]討論的數(shù)據(jù),使用此數(shù)據(jù)有效地證明了新提出的GLET-GTII 分布的優(yōu)越性,這組數(shù)據(jù)包括3 月1 日至3 月30 日期間的歐洲國家因新冠肺炎而每日死亡人數(shù)的觀測結(jié)果,該數(shù)據(jù)在https：//covid19.who.int/中報(bào)告,由30 個(gè)單位組成的數(shù)據(jù)集如下：6、18、29、28、47、55、40、150、129、184、236、237、336、219、612、434、648、706、838、1 129、1 421、118、1 393、1 540、1 941、2 175、2 278、2 824、2 803 和2 667.歐洲國家新冠肺炎數(shù)據(jù)集的描述性度量結(jié)果見表2.

表2 歐洲國家新冠肺炎數(shù)據(jù)集的描述性度量Tab.2 Descriptive degree of covid-19 pneumonia data sets for European countries

通過表2 可知, 歐洲國家在3 月1 日至3 月30 日期間平均每天因新冠肺炎死亡的人數(shù)約841 人,其中最少1 天的死亡人數(shù)是6 人,最多一天的死亡人數(shù)達(dá)到2 818 人.該數(shù)據(jù)的偏度為0.92>0,呈現(xiàn)左偏狀態(tài).標(biāo)準(zhǔn)差為938.69,說明這組數(shù)據(jù)的離散程度較大.圖4 分別繪制了這組數(shù)據(jù)的箱線圖、QQ 圖以及經(jīng)驗(yàn)分布圖,其中箱線圖顯示了這組數(shù)據(jù)集的匯總信息,陰影部分分別展現(xiàn)了25%、50%和75%的分位數(shù),此圖較好地指出了這組數(shù)據(jù)中的值是如何分布的.

圖4 歐洲國家新冠肺炎數(shù)據(jù)集的基本統(tǒng)計(jì)圖Fig.4 Basic statistical graph of the covid-19 data set for European countries

表3 展示了這組真實(shí)數(shù)據(jù)集在六種不同的分布下進(jìn)行參數(shù)估計(jì)的結(jié)果,為了檢查數(shù)據(jù)是否符合MOAPGT-II分布,采用了表4 中的五種評價(jià)準(zhǔn)則.從結(jié)果可以看出,MOAPGT-II 分布擬合當(dāng)前數(shù)據(jù)與選取的五個(gè)測試分布相比,發(fā)現(xiàn)MOAPGT-II 分布的負(fù)對數(shù)似然值為231.90,AIC 為471.79,CAIC為472.75,HQIC 為473.58,AIC、HQIC、CAIC 和負(fù)對數(shù)似然值均小于其他分布,這表明新提出的MOAPGT-II 分布比ETGT-II、GIWD、Frechet、GT-II 以及EB-XII 分布更適合這個(gè)數(shù)據(jù),說明這個(gè)數(shù)據(jù)與MOAPGT-II分布幾乎完美匹配,MOAPGT-II 分布在所有競爭分布中提供了最佳擬合.

表3 歐洲國家新冠肺炎數(shù)據(jù)集在不同分布情況的參數(shù)估計(jì)和標(biāo)準(zhǔn)誤差Tab.3 Parameter estimates and standard errors for different distributions of covid-19 in European national datasets

表4 對于歐洲國家新冠肺炎數(shù)據(jù)集的擬合分布評價(jià)準(zhǔn)則值Tab.4 Evaluation criteria values for the fitted distribution of the European covid-19 data set

4 小結(jié)

本文提出了一種新的MOAPGT-II 分布及其基本性質(zhì).通過使用極大似然估計(jì)方法對該分布的參數(shù)進(jìn)行估計(jì),利用蒙特卡洛模擬分析方法得到的結(jié)果為：隨著樣本量的增大,其平均偏差和均方誤差均會減小并接近0.最后通過一組真實(shí)的歐洲國家因新冠肺炎死亡人數(shù)的案列,將MOAPGT-II 分布與II 型耿貝爾分布、廣義逆Weibull 分布、指數(shù)變換的Gumbel Type-II 分布、Frechet 分布、指數(shù)Burr XII 分布進(jìn)行比較,基于經(jīng)典的擬合優(yōu)度指標(biāo),如AIC、CAIC 以及HQIC,發(fā)現(xiàn)MOAPGT-II 分布擬合該數(shù)據(jù)集的能力要優(yōu)于這些分布.需要進(jìn)一步研究的問題是將對數(shù)法、百分位法[12]、L 矩估計(jì)[13]、矩估計(jì)、EM算法[14-16]以及貝葉斯估計(jì)[17-18]等參數(shù)估計(jì)方法應(yīng)用到本文所提出的分布上,且考慮以后可以使用boostrap 法進(jìn)行區(qū)間估計(jì).