馬永傳，武 東，朱雅敏

(1.皖西學院 金融與數學學院，安徽 六安 237009；2.安徽農業大學 理學院，安徽 合肥 230036；3.華中科技大學 數學與統計學院，湖北 武漢 430074)

?

廣義逐次定數截尾下Pareto分布的統計分析

馬永傳1，武 東2，朱雅敏3

(1.皖西學院 金融與數學學院，安徽 六安 237009；2.安徽農業大學 理學院，安徽 合肥 230036；3.華中科技大學 數學與統計學院，湖北 武漢 430074)

基于廣義逐次定數截尾樣本,對Pareto分布壽命產品進行貝葉斯統計分析,利用吉布斯抽樣給出該模型的近似Bayes估計,最后通過模擬例子表明Bayes估計是有效的.

廣義逐次定數截尾;Pareto分布;Bayes估計;吉布斯抽樣

大量文獻研究了經典的定時截尾和定數截尾及其統計分析[1]．Balakrishnan和Aggarwala(2000)將經典的定時截尾和定數截尾推廣到廣義逐次定時截尾和廣義逐次定數截尾[2]．廣義逐次定數截尾試驗在可靠性分析中經常會遇見，擁有巨大價值，可廣泛應用于醫學、電子技術、機械和航空等領域．在Pareto分布產品的可靠性研究中，王亮等[3]提出了逐步首失效樣本，并就Pareto分布進行了Bayes分析，侯華蕾等[4]考慮雙邊定數截尾壽命試驗，該試驗的優點就是允許左截尾和右截尾，主要利用Bayes方法對雙邊定數截尾下Pareto分布參數進行了估計．衛超等[5]研究了逐步增加II型截尾與定時截尾的融合稱為逐步II型混合截尾，也對Pareto分布進行了Bayes分析．李鳳[6-8]解決了逐步增加II型截尾下Pareto分布參數的逆矩估計、區間估計和Bayes估計．李瓊等[9]也研究了逐步增加II型截尾下Pareto分布的Bayes估計，討論參數先驗兩種方法：共軛先驗與無信息先驗法，并利用吉布斯抽樣法解決Bayes計算問題．但以上均未對廣義逐次定數截尾下Pareto分布的參數估計展開研究，本文就此問題進行討論，主要分為3個部分：①給出廣義逐次定數截尾下的模型進行描述 ；②解決廣義逐次定數截尾下Pareto分布參數的最大似然估計；③討論廣義逐次定數截尾下Pareto分布參數的Bayes估計，結合吉布斯抽樣法給出算法步驟；④利用Monte Carlo方法產生Pareto分布廣義逐次定數截尾樣本，并對Bayes估計與最大似然估計進行精度分析，結果表明Bayes估計有效而實用．

1 模型描述

廣義逐次定數截尾試驗是在定數截尾基礎上的拓展，其模型描述和安排如下：首先依據隨機原則選擇n個測試產品進行試驗，觀察到第一個測試產品失效時刻xr+1:m:n時，已經有r個測試產品失效，然后從所剩的n-r-1個測試產品里隨機移走Rr+1個產品，將余下的n-r-Rr+1-1個測試產品進行試驗；觀測到第二個測試產品失效時刻xr+1:m:n時，從所剩的n-r-Rr+1-2個測試產品里隨機移走Rr+2個產品，將余下的n-r-Rr+1-Rr+2-2個測試產品進行試驗；一直這樣進行下去，直至觀察到預定的m-r個測試產品失效時刻xm:m:n時結束試驗．最后將所剩的Rm=n-r-Rr+1-Rr+2-…-Rm-1-(m-r-1)個測試產品移走．容易知道m-r個產品的失效時間為xr+1:m:n≤xr+2:m:n≤…≤xm:m:n.依據文獻[2]知道似然函數如下：

(1)

2 最大似然估計

假設1 假設測試產品的壽命服從Pareto分布P(μ,α)，其概率密度和分布函數分別為

(2)

和

(3)

其中μ為分布的門限參數，α為分布的刻度參數．

設xr+1:m:n≤xr+2:m:n≤…≤xm:m:n是來自Pareto分布樣本容量為n的廣義逐次定數截尾樣本，Rr+1,Rr+2,…,Rm-1,Rm是試驗中依次被移開的產品個數．將(2)和(3)式代入(1)式，可得廣義逐次定數截尾下Pareto分布的似然函數為

(4)

這里data={xr+i:m:n|i=1,2,…,m-r}．

于是，基于Pareto分布的廣義逐次定數截尾樣本的對數似然函數為

(5)

通過求解上述似然方程組，可得參數μ和α的最大似然估計分別為

和

3 Bayes估計

從而得到參數μ的滿條件后驗密度為

而參數α的滿條件后驗密度為

下面采用吉布斯抽樣(Gibbs Sampler)，對參數進行模擬和估計．它們隨機數的產生都可以采用舍選抽樣法．下面給出了廣義逐次定數截尾下Pareto分布的Bayes估計的步驟：

Step1：選取參數(μ.α)的起始點(μ(0),α(0))，一般采用先驗分布產生；

Step2：依據參數α的滿條件后驗密度π(α|μ(i-1))，并從中抽取隨機數α(i)；

Step3：依據參數μ的滿條件后驗密度π(μ|α0(i))，并從中抽取隨機數μ(i)；

Step4：令i=i+1，然后返回Step2，直至達到給定的迭代次數便結束迭代．

則α0(i)，μ(i)，i=1,2,…,N為參數(μ，α)的一個吉布斯迭代樣本，通常將前N0(N>N0)個吉布斯迭代樣本去掉，用剩余的N-N0個迭代樣本的期望值作為參數(μ，α)的Bayes估計，即(μ,α)的Bayes估計分別為:

4 Monte Carlo模擬

以上對廣義逐次定數截尾下Pareto分布進行了Bayes分析，利用Monte Carlo方法可得到該模型的隨機數，具體步驟如下：

(1)從標準指數分布產生r+1個相互獨立的隨機數W1,W2,…,Wr+1；

(3)從標準指數分布產生m-r-1個相互獨立的隨機數Zr+2,Zr+3,…,Zm；

(4)令

i=2,3,…，m-r,則得到的廣義逐次定數截尾下標準指數分布樣本；其中n為試驗的產品總個數，Rr+1,Rr+2,…,Rm-1為試驗中依次被撤離的測試產品個數．

為了考察最大似然估計(MLE)與Bayes估計(Bayes)的精度，現取Pareto分布參數為α=1.5,σ=300，從該分布產生m-r個廣義逐次定數截尾樣本xr+1,xr+2,…,xm．每種方案產生100組模擬樣本并分別計算估計的相對偏差與均方誤差，這里Bayes估計迭代步數N=11000,舍棄迭代步數為N0=1000，迭代初始值為α=3,σ=200，具體試驗安排與估計結果如表1所示．

表1 廣義逐次定數截尾下Pareto分布的試驗方案與估計結果

從表1可以看出，該模型Bayes估計的精度較高，說明估計更加有效．

5 結語

利用吉布斯抽樣法獲得了廣義逐次定數截尾下Pareto分布場合的一種近似Bayes估計．關于先驗分布的選取，Pareto分布的門限參數和形狀參數的先驗取為無信息先驗分布，而吉布斯抽樣迭代過程的隨機數產生采用了舍選抽樣法．在基于吉布斯抽樣法的Bayes計算過程，需要產生大量隨機數，這樣會造成計算效率低下,可考慮采用基于Lindley公式的積分近似法計算后驗期望．

[1]茆詩松,湯銀才,王玲玲.可靠性統計[M].北京:高等教育出版社,2008.

[2]BALAKRISHNAN N,AGGARWALA R.Progressive censoring:Theory,methods and applications[M].Boston: Birkhguser,2000.

[3]王亮,師義民,孫天宇.兩參數Pareto分布逐步首失效樣本的Bayes估計[J].系統工程理論與實踐,2012,32(11): 2498-2503.

[4]侯華蕾,師義民,李豪亮.雙邊定數截尾下Pareto分布的可靠性分析[J].數理統計與管理,2009,28(5):826-830.

[5]衛超,師義民.逐步II型混合截尾下Pareto分布的統計分析[J].火力與指揮控制,2014,39(4):27-30.

[6]李鳳.逐次定數截尾下Pareto分布參數的逆矩估計[J].統計與決策,2010(24):156-157.

[7]李鳳.逐步增加II型截尾下Pareto分布參數的區間估計[J]. 科學技術與工程,2011,11(11):2554-2557.

[8]李鳳,師義民.逐步增加II型截尾下Pareto分布的Bayes估計[J].數學的實踐與認識,42(13):137-142.

[9]李瓊,武東.逐步增加定數截尾下Pareto分布的貝葉斯分析[J].上海第二工業大學學報,2012,29(2):135-138.

(責任編輯:陳衍峰)

2016-05-06

國家自然科學基金項目“面板數據的最優設計”(11401056)

馬永傳,安徽六安人,副教授；武東,安徽六安人,副教授；朱雅敏,女,安徽安慶人．

O213.2

A

1008-7974(2016)05-0022-03

10.13877/j.cnki.cn22-1284.2016.10.007