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

王丙參，魏艷華

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

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

王丙參，魏艷華

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

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

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

0 引言

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

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

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

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

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

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

設(shè)隨機(jī)變量Xj表示T在區(qū)間[τj-1,τj)上的截?cái)喾植迹琷=1,2,…,k，則Xj的條件密度為：

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

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

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

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

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

②計(jì)算接受候選樣本α′的概率

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

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

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

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

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

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

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

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

α的完全條件分布

λ的完全條件分布

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

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

②計(jì)算

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

①?gòu)腘(λk,σ2λ)抽取樣本λ′，若λ′≤0，則重新抽樣。

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

2 蒙特卡羅模擬

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

例1：設(shè)壽命T服從廣義指數(shù)分布，形狀參數(shù)α=2，尺度參數(shù)λ=1，根據(jù)重點(diǎn)抽樣原則，在密度函數(shù)取值大的區(qū)域劃分更細(xì)點(diǎn)，不妨取檢測(cè)時(shí)刻τ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ù)（見(jiàn)表1），并給出模擬過(guò)程中目標(biāo)參數(shù)的軌跡圖（左）、直方圖（中）和自相關(guān)系數(shù)圖（右），本文只給出第1組數(shù)據(jù)的結(jié)果（見(jiàn)圖1），其他類似。

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

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

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

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

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

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

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

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

(2.0336,0.0825)，(1.0161,0.0162)

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

3 結(jié)論

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

[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ì)[J].系統(tǒng)工程理論與實(shí)踐,2011,31(9).

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

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

[6]Givens G H,Hoeting J A.計(jì)算統(tǒng)計(jì)[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ù)場(chǎng)合逆威布爾分布參數(shù)貝葉斯估計(jì)的混合Gibbs算法[J].四川大學(xué)學(xué)報(bào):自然科學(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é)院中青年教師科研資助項(xiàng)目（TSA1506）

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

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