程貝麗,周菊玲
(新疆師范大學 數(shù)學科學學院,新疆 烏魯木齊 830017)
對數(shù)伽瑪分布是一種重要的壽命分布,廣泛應用在生物學、金融及新型農(nóng)村合作醫(yī)療住院費用等可靠性領域[1-4].由于對數(shù)伽瑪分布具有向右偏斜的特點,所以在非壽險精算中也得到了一定的應用.Zhang等[1]從生物學角度研究了對數(shù)伽瑪分布的模型;Amini等[5]研究了對數(shù)伽瑪分布族的性質;王成元等[6-8]研究了不同損失函數(shù)下對數(shù)伽瑪分布尺度參數(shù)的貝葉斯估計,并對不同的貝葉斯估計做出比較;金秀巖[9-10]分別探討了對數(shù)伽瑪分布和負對數(shù)伽瑪分布條件概率的封閉性、Pearson-χ2距離和漸進性;熊福生[11]研究了對數(shù)伽瑪分布和負對數(shù)伽瑪分布的再生性,得到了一些與再生性、準再生性相關的結果.在近幾年的研究中,變點問題一直是統(tǒng)計學中的熱點問題.變點問題一般的研究方法有最小二乘法、貝葉斯方法、極大似然方法、累次計數(shù)法和非參數(shù)法.陳希孺[12]結合概率變點問題,介紹了極大似然音樂方法、累計次數(shù)法和貝葉斯方法;黃月蘭[13]分別用極大似然方法、貝葉斯方法和最小二乘法討論了單參數(shù)指數(shù)分布的變點問題.沙雪云等[14]用極大似然方法和貝葉斯方法討論了單參數(shù)Parato分布的變點問題;陳杰[15]分別在非貝葉斯場景下和貝葉斯場景下對變點進行檢測.但是,用不同的方法研究對數(shù)伽瑪分布變點問題的文章還比較少,本文用極大似然估計和貝葉斯估計研究單變點對數(shù)伽瑪分布的參數(shù)估計問題.
下面分別用極大似然估計方法和貝葉斯估計方法對單變點對數(shù)伽瑪分布的參數(shù)進行估計.
設隨機變量X1,X2,…,Xn為來自對數(shù)伽瑪分布模型中的一個樣本,則對其n個獨立觀測值x1,x2,…,xn得到的似然函數(shù)為
(1)
令γ=(θ1,θ2,k),則單變點對數(shù)伽瑪分布的似然函數(shù)為
(2)
固定k,對θ1,θ2求極大,得
(3)
在統(tǒng)計推斷中,樣本的狀態(tài)固然重要,但歷史的經(jīng)驗也同樣舉足輕重.故為了提高統(tǒng)計推斷的質量,應注意先驗信息對整個統(tǒng)計推斷的影響,沒有先驗分布或取不恰當?shù)南闰灧植级伎赡苁菇y(tǒng)計推斷的結果產(chǎn)生誤差.
首先確定各參數(shù)的先驗分布:
1)變點k的先驗分布為無信息先驗分布:π(k)∝1;
2)尺度參數(shù)θ1,θ2的先驗分布分別取無信息先驗分布和共軛先驗分布.
討論尺度參數(shù)θ1,θ2的先驗分布取無信息分布時參數(shù)的貝葉斯估計,假定θ1,θ2,k之間相互獨立.
引理1 設隨機變量X服從對數(shù)伽瑪分布,參數(shù)θi的先驗分布π(θi)為無信息先驗分布,X1,X2,…,Xn為服從對數(shù)伽瑪分布的一個樣本,記X=(X1,X2,…,Xn),T=(t1,t2),則θi的后驗分布服從Ga(nα,T).
證由式(2)可知,x1,x2,…,xn是樣本X1,X2,…,Xn的n個獨立觀測值的似然函數(shù),又已知參數(shù)θi的先驗分布為無信息先驗分布,即π(θi)∝1.則參數(shù)θi的后驗概率密度函數(shù)為
于是θi|x~Ga(nα,ti),故θi|X~Ga(nα,T).
則當參數(shù)θi的先驗分布為無信息先驗分布時,對數(shù)伽瑪分布單變點模型的似然函數(shù)為
(4)
得到各參數(shù)的滿條件分布如下:
討論尺度參數(shù)θ1,θ2的共軛先驗分布為伽瑪分布時參數(shù)的貝葉斯估計,假定θ1,θ2,k之間相互獨立.
引理2[3]設隨機變量X服從對數(shù)伽瑪分布,參數(shù)θi的先驗分布π(θi)為伽瑪分布Ga(ai,bi),X1,X2,…,Xn是服從對數(shù)伽瑪分布的一個樣本,記X=(X1,X2,…,Xn),T=(t1,t2),則θi的后驗分布服從Ga(nα+ai,T+bi).
則參數(shù)θi的后驗概率密度函數(shù)為
于是θi|x~Ga(nα+ai,bi+ti),故θi|X~Ga(nα+ai,bi+T).
則當θi的先驗分布為伽瑪分布時,對數(shù)伽瑪分布單變點模型的似然函數(shù)為
(5)
得到各參數(shù)的滿條件分布如下:
尺度參數(shù)θ1,θ2滿條件分布的形式較為簡單,可直接用Gibbs抽樣,而變點k的滿條件分布較為復雜,不再適合Gibbs抽樣方法,則利用MH算法進行抽樣.
MCMC算法可分為如下幾個具體步驟.
為了確保抽樣的收斂性,在模擬過程中,先進行10 000次的預迭代,再進行20 000次的Gibbs迭代,即從10 001次開始到30 000次結束.假設Gibbs樣本的容量為M,如果第B次以后逐漸收斂,則把后M-B個樣本的均值作為各未知參數(shù)的估計值.
根據(jù)對數(shù)伽瑪分布的概率密度函數(shù)和極大似然估計式(3),利用MATLAB軟件對未知參數(shù)θ1,θ2,k進行估計,結果如表1所示.
表1 參數(shù)θ1,θ2,k的極大似然估計
模擬結果顯示,極大似然方法能夠估計未知參數(shù)的值,但參數(shù)θ1,θ2,k的極大似然估計值與真值相差較大,精度不高.
取超參數(shù)ai=bi=1,i=1,2,取B=10 000,M=30 000,根據(jù)各參數(shù)的滿條件分布進行MCMC模擬,將20 000個樣本的均值作為各未知參數(shù)的估計值.
1)當先驗分布取無信息先驗分布時,參數(shù)θ1,θ2,k的參數(shù)估計如表2所示.
表2 無信息先驗分布情況下參數(shù)θ1,θ2,k的參數(shù)估計
2)當先驗分布取共軛先驗分布時,參數(shù)θ1,θ2,k的參數(shù)估計如表3所示.
表3 共軛先驗分布情況下參數(shù)θ1,θ2,k的參數(shù)估計
圖1、圖3是變點k抽樣迭代圖,是判斷Gibbs樣本是否收斂的重要依據(jù).用MCMC方法模擬時,會同時產(chǎn)生兩條迭代鏈,若這兩條迭代鏈逐漸穩(wěn)定并趨于重合,則說明Gibbs抽樣是收斂的,變點k的兩條迭代鏈如圖2、圖4所示.
圖1 變點k的抽樣迭代圖 圖2 變點k的兩條迭代鏈
圖3 變點k的抽樣迭代圖 圖4 變點k的兩條迭代鏈
下面對模擬的結果進行分析,從表1、表2中可以看到θ1,θ2,k這3個參數(shù)的均值與真值都較為接近,且MC誤差都小于2%,即各參數(shù)的貝葉斯估計精度較高;各參數(shù)的標準差也比較小,說明數(shù)據(jù)比較穩(wěn)定;各參數(shù)置信水平0.95的區(qū)間[2.5%分位數(shù),97.%分位數(shù)]的長度也較短,即得到的區(qū)間估計效果較好.由圖1和圖3可以看出Gibbs抽樣基本在變點附近波動,呈現(xiàn)出一定的規(guī)律性.其次,觀察圖2和圖4,可以發(fā)現(xiàn)產(chǎn)生的兩條馬爾科夫鏈穩(wěn)定且趨于重合,說明收斂性較好.綜上所述,對數(shù)伽瑪分布變點問題可以用MCMC方法實現(xiàn),且在不同先驗分布下,估計值的精度不同.
本文主要研究對數(shù)伽瑪分布尺度參數(shù)存在變點時的參數(shù)估計問題.隨機模擬的結果表明:兩種方法都能夠有效地對未知參數(shù)進行估計.在相同條件下,貝葉斯估計所需要的時間遠遠大于極大似然估計,但是相比于極大似然估計,貝葉斯估計的值更穩(wěn)定,也更接近真實值.通過比較可知,共軛先驗分布下的貝葉斯估計值優(yōu)于無信息先驗下的貝葉斯估計值.這說明了選取恰當?shù)南闰灧植伎梢蕴岣呓y(tǒng)計推斷的質量.因此,基于本文中單變點對數(shù)伽瑪分布的模型,先驗分布取共軛先驗分布時,未知參數(shù)的貝葉斯估計是最優(yōu)的.