王丙參，魏艷華

（天水師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，甘肅 天水 741001）

不完全數(shù)據(jù)場合廣義指數(shù)分布參數(shù)Bayes估計的混合Gibbs算法

王丙參，魏艷華

（天水師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，甘肅 天水 741001）

文章針對在兩種不完全數(shù)據(jù)場合，給出了廣義指數(shù)分布參數(shù)貝葉斯估計的混合Gibbs算法，Monte-Carlo模擬結(jié)果顯示：混合Gibbs算法簡單、可行、精度較高、適應(yīng)范圍廣。分組數(shù)據(jù)的分組方式對模擬結(jié)果影響較大。

分組數(shù)據(jù);定數(shù)截尾樣本;廣義指數(shù)分布;Bayes估計;混合Gibbs算法

0 引言

廣義指數(shù)分布是Gupta和Kundu于1999年提出的，作為指數(shù)分布的一種推廣，它在可靠性理論中有廣泛應(yīng)用[1-3]。對于實際問題，人們根據(jù)經(jīng)驗往往知道一定的先驗信息，所以目標(biāo)量的貝葉斯估計會比相對經(jīng)典估計的效果好。一般情況下，目標(biāo)量的貝葉斯估計可轉(zhuǎn)化為關(guān)于被積函數(shù)π(x)（即后驗分布）的積分，但在不完全數(shù)據(jù)場合，π(x)可能很復(fù)雜，無法用單一的Gibbs抽樣獲得其貝葉斯估計。將Gibbs算法和自適應(yīng)算法方法結(jié)合通常推導(dǎo)繁瑣，模擬效果一般且受后驗分布影響大[4,5]。因此，希望找到一種簡單、穩(wěn)定、效果好的抽樣算法。因為Gibbs算法可將高維抽樣轉(zhuǎn)化為低維抽樣（常用一維抽樣），而Metropolis算法特別適合低維抽樣，所以將二者結(jié)合，可達(dá)到出人意料的效果[6-8]。本文在不完全數(shù)據(jù)場合將Gibbs算法和Metropolis算法混合（即在Gibbs抽樣中使用Metropolis算法抽樣）求廣義指數(shù)分布參數(shù)的貝葉斯估計。

1 廣義指數(shù)分布參數(shù)貝葉斯估計的混合Gibbs算法

由于產(chǎn)品壽命受多種因素（如使用頻率、環(huán)境等）影響，故可用一個隨機(jī)變量表示。假如某種產(chǎn)品壽命T服從廣義指數(shù)分布,即分布函數(shù)（cdf）和密度函數(shù)（pdf）分別為：

1．1 分組數(shù)據(jù)場合

在初始時刻τ0=0同時隨機(jī)投入n個產(chǎn)品試驗，由于各種（如成本、條件限制等）原因，我們只能定期檢測，獲得分組數(shù)據(jù)：假定檢測時刻τ0,τ1,…,τk滿足0=τ0＜τ1＜… ，＜τk-1＜τk=∞r(nóng)j為時間區(qū)間[τj-1,τj)上的失效個數(shù)，j=1,2,…,k，失效時間t1,t2,…,tn只知落入?yún)^(qū)間，具體數(shù)值未知，Y為已觀測情況。

對于分組數(shù)據(jù)，由于無法獲得聯(lián)合分布f(t1,…,tn|α,λ,Y)的樣本，所以需要完全條件分布的Gibbs抽樣難于操作。于是，先令t=(t1,t2,…,tn)，其中t1≤t2≤…≤tn為具體的失效時間，再考察(t,α,λ)的后驗分布函數(shù)。顯然，t的完全條件分布：

設(shè)隨機(jī)變量Xj表示T在區(qū)間[τj-1,τj)上的截斷分布，j=1,2,…,k，則Xj的條件密度為：

于是，分別從截斷分布Xj中抽取rj個樣本，j=1,2,...,k，便可得到n個聯(lián)合分布f(t1,…,tn|α,λ,Y)的樣本，不妨記為t1,…,tn。這樣，樣本t1,t2,…,tn的聯(lián)合pdf可表示為：

當(dāng)形狀參數(shù)α與尺度參數(shù)λ未知時，如果取參數(shù)α,λ的聯(lián)合先驗密度 π(α,λ)∝ 1,則其聯(lián)合后驗pdf為：

形狀參數(shù)α與尺度參數(shù)λ的完全條件分布分別為：

顯然，對參數(shù)α,λ直接抽樣比較困難，下面嘗試用混合Gibbs算法進(jìn)行抽樣，具體操作如下：在Gibbs抽樣中，假定初值 (α(0),λ(0))給定，(α(i),λ(i))表示第i次迭代估計值，它可作為第i+1次迭代的初始值，則第i+1次Gibbs迭代可分為3步。

①從α的建議分布抽取樣本α′，若α′≤0，則重新抽樣，其中αk為當(dāng)前狀態(tài)，σα為為建議分布標(biāo)準(zhǔn)差。σα用來控制候選樣本α′的波動，取值要合適[6]。

②計算接受候選樣本α′的概率

③生成隨機(jī)數(shù)u～U(0,1)，令

④令k=k+1，返回第①步。

③生成隨機(jī)數(shù)u～U(0,1)，令

④令k=k+1，返回第①步。

在（2）、（3）步中，利用Metropolis算法可得到一個馬氏鏈，當(dāng)它達(dá)到均衡狀態(tài)后，馬氏鏈的任一個樣本就可作為α(i+1)與λ(i+1)，這便完成一次迭代：令i=i+1，重復(fù)上述（1）—（3）步即可得α,λ的后驗分布樣本并進(jìn)行貝葉斯估計。當(dāng)然，間隔取樣可在一定程度上消減后驗分布樣本的相關(guān)性，從而保證樣本的近似獨立性，間隔大小取決于樣本的延遲自相關(guān)系數(shù)。顯然，混合Gibbs算法很容易推廣到更多參數(shù)情況。由于混合Gibbs算法一般對參數(shù)α,λ的不同先驗分布都較容易抽樣，故本文只給出上述一種先驗分布的情形，其他先驗分布情況類似。

1．2 定數(shù)截尾樣本場合

設(shè)隨機(jī)抽取的n個產(chǎn)品中有r個失效時便停止試驗，得次序失效數(shù)據(jù)t1≤t2≤…≤tr，其聯(lián)合pdf：

同上，仍取α,λ的聯(lián)合先驗分布 π(α,λ)∝ 1,則它們聯(lián)合后驗分布：

α的完全條件分布

λ的完全條件分布

同分組數(shù)據(jù)一樣，這里仍需使用混合Gibbs算法抽樣，具體步驟為：在Gibbs抽樣的第i+1次迭代中，分為以下2步。

（1）利用M算法從 π(α|t1,t2,…,tr,λ(i))抽取樣本α(i+1)。

②計算

（2）利用M算法從 π(λ|t1,t2,…,tr,α(i+1))抽取樣本λ(i+1)。

①從N(λk,σ2λ)抽取樣本λ′，若λ′≤0，則重新抽樣。

求貝葉斯估計的具體做法同上，不再詳述。

2 蒙特卡羅模擬

利用MATLAB（2014a）在不完全數(shù)據(jù)場合對廣義指數(shù)分布的混合Gibbs算法進(jìn)行隨機(jī)模擬。

例1：設(shè)壽命T服從廣義指數(shù)分布，形狀參數(shù)α=2，尺度參數(shù)λ=1，根據(jù)重點抽樣原則，在密度函數(shù)取值大的區(qū)域劃分更細(xì)點，不妨取檢測時刻τ0=0,τ1=0.45,τ2=0.85,τ3=1.15,τ4=1.55,τ5=1.95,τ6=2.95,τ7=4,τ8= ∞ ，由于落在 [τj-1,τj)的樣本量rj不同，j=1,2,…,8，本文利用反函數(shù)法生成廣義指數(shù)分布隨機(jī)數(shù)，樣本容量n=200，依次取得下面5組分組數(shù)據(jù)（見表1），并給出模擬過程中目標(biāo)參數(shù)的軌跡圖（左）、直方圖（中）和自相關(guān)系數(shù)圖（右），本文只給出第1組數(shù)據(jù)的結(jié)果（見圖1），其他類似。

表1 廣義指數(shù)分布五組分組數(shù)據(jù)

圖1 第1組分組數(shù)據(jù)的Bayes估計結(jié)果

對于這5組數(shù)據(jù)，均選取初值α=2.5，λ=1.5，建議分布為正態(tài)分布，其期望為參數(shù)當(dāng)前值，標(biāo)準(zhǔn)差為0.4。為保證馬氏鏈達(dá)到穩(wěn)定狀態(tài)，去除開始的100迭代。同理，在每一次Gibbs抽樣中的M算法采用第101次結(jié)果作為樣本。由自相關(guān)圖（圖1右）可知，在Lag=5時樣本的自相關(guān)系數(shù)很小，所以每隔5次迭代取1個樣本可近似保證樣本的獨立性，共得到α,λ的各1000個樣本，即迭代5000次。這1000個樣本均值就是參數(shù)在平方損失函數(shù)下的貝葉斯估計，利用樣本分位數(shù)可求得參數(shù)的95%置信區(qū)間，見表2。

表2 廣義指數(shù)分布參數(shù)的貝葉斯估計與置信區(qū)間

由表2可知，利用混合Gibbs算法求得廣義指數(shù)分布參數(shù)的α,λ貝葉斯估計，估計均值比較接近真值，置信區(qū)間均覆蓋真實值且相對較短。采用5組估計值的均值作為最終結(jié)果：=1.9930，=0.9995，其均方誤差分別為0.0761，0.0569，即均方誤差較小，模擬結(jié)果相對穩(wěn)定。進(jìn)一步模擬可知，檢測時間取值不同會對模擬結(jié)果有較大影響，因此，在采集分組數(shù)據(jù)時也要深入思考，盡量選取合理分組。樣本容量對模擬結(jié)果也影響顯著，特別對于不完全數(shù)據(jù)場合，如果樣本量少，則誤差可能較大。

例2：設(shè)壽命T服從參數(shù) (α,λ)為(2,1)和(1,1)的廣義指數(shù)分布，在20%定數(shù)截尾情況下蒙特卡羅模擬，每次模擬均利用反函數(shù)法產(chǎn)生200個廣義指數(shù)分布隨機(jī)數(shù)。同上，剔除前100次迭代，每隔5次迭代取1個樣本，共獲得目標(biāo)參數(shù)1000個樣本。每次Gibbs抽樣中的M算法取第101次迭代結(jié)果作為參數(shù)的樣本。部分模擬結(jié)果見表3，目標(biāo)參數(shù)的軌跡圖、直方圖和自相關(guān)系數(shù)圖省略。

表3 1次模擬參數(shù)的貝葉斯估計與95%置信區(qū)間

對于 (α,λ)=(2,1)，下面進(jìn)行100次蒙特卡羅模擬，目標(biāo)參數(shù)(α,λ)的期望與標(biāo)準(zhǔn)差分別為：

(2.0336,0.0825)，(1.0161,0.0162)

這表明在定數(shù)截尾樣本場合，利用混合Gibbs算法求得廣義指數(shù)分布的Bayes估計精度較高。

3 結(jié)論

對于不完全數(shù)據(jù)場合的廣義指數(shù)分布參數(shù)的貝葉斯估計，本文采用的混合Gibbs算法理論簡單，約束少，對不同的先驗分布都容易實現(xiàn)，且精度較高、穩(wěn)定，非常適合具有復(fù)雜后驗分布的情形。模擬結(jié)果顯示，分組數(shù)據(jù)的不同分組會對估計結(jié)果有較大影響，因此分組要慎重。利用模擬隨機(jī)數(shù)代替樣本的具體觀測值是處理不完全數(shù)據(jù)場合（如缺失數(shù)據(jù)場合、定時截尾等）參數(shù)估計的常用手段。

[1]Gupta R D,Kundu D.Generalized Exponential Distribution[J].Austral&New Zealand J Statist,1999,41(2).

[2]Gupta R D,Kundu D.Generalized Exponential Distribution:Existing Results and Some Recent Developments[J].Journal of Statistical Planning&Inferencem 2007,137(11).

[3]郡偉安,師義民,孫天宇等.熵?fù)p失函數(shù)下廣義指數(shù)分布的可靠性估計[J].系統(tǒng)工程理論與實踐,2011,31(9).

[4]湯銀才,侯道燕.三參數(shù)Weibull分布參數(shù)的Bayes估計[J].系統(tǒng)科學(xué)與數(shù)學(xué),2009,29(1).

[5]劉飛,王祖堯,竇毅芳等.基于Gibbs抽樣算法的三參數(shù)威布爾分布Bayes估計[J].機(jī)械強(qiáng)度,2007,29(3).

[6]Givens G H,Hoeting J A.計算統(tǒng)計[M].王兆軍,劉民千等譯.北京:人民郵電出版社,2009.

[7]Altaleb A,Chauveau D.Bayesian Analysis of the Logit Model and Comparison of Two Metropolis-hastings Strategies[J].Computational Statistics&Data Analysis,2002,39(2).

[8]魏艷華,王丙參,孫永輝.分組數(shù)據(jù)場合逆威布爾分布參數(shù)貝葉斯估計的混合Gibbs算法[J].四川大學(xué)學(xué)報:自然科學(xué)版,2014,51(1).

Mixed Gibbs Algorithm for Bayesian Estimation of the Parameters of Generalized Exponential Distribution in Incomplete Data Occasion

Wang Bingcan，Wei Yanhua

(School of Mathematics and Statistics,Tianshui Normal University,Tianshui 741001,China)

This paper gives a mixed Gibbs algorithm of Bayesian estimation of the parameters of generalized exponential distribution in the two occasions of incomplete data.The Monte-Carlo simulation result shows that the mixed Gibbs algorithm is simple,feasible and accurate with a wide range of adaptability,and that the grouping mode of packet data has a great influence on the simulation results.

grouped data;fixed censored samples;generalized exponential distribution;Bayesian estimation;mixed Gibbs algorithm

O213

A

1002-6487（2017）20-0023-03

天水師范學(xué)院中青年教師科研資助項目（TSA1506）

王丙參（1983—），男，河南南陽人，碩士，講師，研究方向：統(tǒng)計計算。

(責(zé)任編輯/亦 民)