逐步增加II型截尾下Pareto分布的贝叶斯分析

李 琼，武 东

（1. 徽商职业学院电子信息系，合肥 230022；2. 安徽农业大学理学院，合肥 230036）

逐步增加II型截尾下Pareto分布的贝叶斯分析

李 琼1，武 东2

（1. 徽商职业学院电子信息系，合肥 230022；2. 安徽农业大学理学院，合肥 230036）

对Pareto分布场合逐步增加II型截尾样本进行了贝叶斯分析，利用马尔可夫链蒙特卡罗方法给出了参数的贝叶斯估计。最后，通过蒙特卡罗模拟和应用实例表明该贝叶斯估计是有效的。

逐步增加II型截尾；Pareto分布；先验分布；贝叶斯估计；Gibbs抽样

0 引言

Pareto分布的失效率函数具有递减性，现常用来研究电子、机械、航天等产品的寿命分布和病人的生存时间等。近年来，在可靠性与生存分析理论中关于Pareto分布统计推断的文献较多。文献[1]和文献[2]分别给出了逐步增加II型截尾下Pareto分布的逆矩估计和区间估计，文献[3]和文献[4]对移除部件数分别服从均匀分布和二项分布两种情形逐步增加II型截尾下Pareto分布进行了统计分析及贝叶斯推断，文献[5]在非对称损失下给出了一般逐步增加截尾下Pareto分布的区间估计。但上述均未讨论采用Gibbs抽样[6]解决逐步增加II型截尾下Pareto分布的贝叶斯统计分析。本文主要分为：第1节，给出了Pareto分布逐步增加II型截尾试验的基本假设与统计模型；第2节，讨论了Pareto分布逐步增加II型截尾试验的两参数的先验选取和利用马尔可夫链蒙特卡罗方法获得的参数的贝叶斯估计；第3节，利用蒙特卡罗方法对逐步增加II型截尾下Pareto分布模型进行了随机模拟，统计表明该贝叶斯估计是有效的；第4节，对HIV/AIDS病人的生存寿命数据进行了逐步增加II型截尾，并给出了基于无信息先验的贝叶斯分析。

1 模型

为了进行逐步增加II型截尾试验，现从一批产品中随机抽取n个产品，按下述方案进行操作：当观察到第一个产品失效时刻X1时，在剩下的n−1个未失效的产品中任意抽取r1个移离试验，将剩下n−1−r1个未失效的产品留下继续进行试验；当观察到第二个产品失效时刻X2时，在剩下的n−1−r1个未失效的产品中任意抽取r2个移离试验，将剩下n−1−r1−r2个未失效的产品留下继续进行试验；如此下去，直到观察到m个产品失效时刻Xm时终止试验，则此时余下rm=n−r1−r2−…rm−1−m 个产品被移离试验。设x1≤x2≤…≤xm是容量为n的逐步增加II型截尾样本，则由文献[7]可得似然函数为

分别为相应的密度函数和分布函数。

若产品服从两参数(,)α σ的Pareto分布，记为Pareto(,)α σ，其密度函数和分布函数分别为和

其中α为形状参数, σ为刻度参数。

设x1≤x2≤…≤xm是来自Pareto分布容量为n的逐步增加II型截尾样本, r1, r2,…,rm是逐步增加II型截尾寿命试验中依次被移离的产品个数。将(2)式和(3)式代入(1)式，得到逐步增加II型截尾样本下Pareto分布的似然函数为

其中data={x1, x2,…,xm}。

容易获得σ和α的最大似然估计分别为

2 贝叶斯估计

为了获得参数的贝叶斯估计，下面利用马尔可夫链蒙特卡罗方法给出逐步增加II型截尾下两参数Pareto分布的一种近似贝叶斯估计的方法。

对产品寿命服从两参数(,)α σ的Pareto分布进行逐步增加II型截尾寿命试验，得到样本数据为x1≤x2≤…≤xm, r1, r2,…,rm是试验中相应被移离的产品数。参数(α ,σ)的联合先验分布取幂伽玛分布[3-4]PG(,,,) ν λµβ, 即

其中ν, λ,µ,β>0,0<βλ<µ。此先验实际上是(α ,σ)的共轭先验分布[3-4]。

根据贝叶斯定理，由似然函数(4)式和先验密度(5)式，得到(α, σ)的联合后验分布为

该分布为伽玛分布Ga(m+ν+1,A), 其随机数可由软件直接生成。

由(6)式可得σ的边际后验分布为

该分布的随机数可按舍选抽样法抽取。

无信息先验（non-informative prior, NIP）能反映参数所有值不受干扰影响时的先验，因此基于无信息先验的贝叶斯推断被认为仅在理论上有价值。然而，当很难获得经典解法的情况下，用无信息先验来进行贝叶斯推断是十分重要的，因为它被视为获得经典方法的一种数学方法。对于Pareto分布，Nigm和Hamdy[8]对两参数(,)α σ使用的无信息先验为

当参数取ν=−1,λ=0,µ=1,β→∞时, (5)式退化为(7)式的形式，即此无信息分布可以看成是共轭先分布的特殊形式。

下面利用马尔可夫链蒙特卡罗方法给出逐步增加II型截尾下Pareto分布的贝叶斯估计，具体可按以下步骤进行：

第一步，假设给出起始点α(0),σ(0);

第二步，从边际后验分布π(α| σ(k−1),data)抽取随机数α(k)；

第三步，从边际后验分布π(σ| α(k),data)抽取随机数σ(k)；

第四步，令k:=k+1，并返回到第二步，直到达到给定终止的迭代次数。

则α(k),σ(k), k=1,2,…,N 为参数(α, σ)的一个Gibbs迭代样本，将前N(N>N)个Gibbs迭代样本舍

00弃，用舍弃后的N−N0个Gibbs迭代样本的样本均值作为(α, σ)的贝叶斯估计，即α, σ的贝叶斯估计分别为

3 随机模拟

利用Gibbs抽样可以获得逐步增加II型截尾下Pareto分布的贝叶斯估计。现在利用Monte Carlo方法进行模拟比较，选取Pareto分布参数为α=1.5,σ=300, 从该Pareto分布产生m个逐步增加II型截尾样本X1,X2,…,Xm, 其具体操作方法可参见文献[9]。

下面利用模拟方法考察无先验信息贝叶斯估计（Bayes）与最大似然估计（MLE）的精度。每种方案产生100组模拟样本分别计算估计量的相对偏差和均方误差，其中贝叶斯估计中迭代次数N=1100，舍弃迭代数取N0=100，先验分布取µ=1,β=5 000,ν=3,λ=0,起始点α=2,σ=200, 试验方案与模拟结果见表1。

从表1可以看出，以上两种估计都达到了较高的精度，但贝叶斯估计比最大似然估计的精度相对较高。

表1 逐步增加II型截尾下Pareto分布的试验方案与估计结果Tab. 1 Tests scheme and estimation results for Pareto distribution under progressive Type-II censoring

4 应用例子

下面利用实际数据研究贝叶斯估计在逐步增加II型截尾试验中的应用，其描述的是100个HIV/AIDS病人的剩余寿命[4]（单位：月），其中1*表示剩余寿命至少在1个月以上，数据如下。

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1*1*2 2 2 2 2 2*2*2*2*2*3 3 3 3 3 3 3 3 3 3 3*3*4 4 4 4 4*5 5 5 5 5 5 5 6 6 6*7 7 7 7 7 7 7*8 8 8 8 9 9 9 10 10 10 10*11 11 11 12 12 12*12*13 14 15 15 19*22 24*30 31 32 34 35 36 43 53 54 56*57 58 60*60*

使用逐步增加II型截尾试验方案S5：n=100,m=12, (r1, r2,…,rm) = (16, 9, 11, 11, 2, 13, 6, 5, 2, 9, 2, 2), 得到逐步增加II型截尾样本为1, 2, 3, 4, 6, 7, 10, 12, 15, 22, 54和58。

文献[4]在平方损失函数下获得逐步增加II型截尾下Pareto分布两参数的贝叶斯估计为σ^=0.979,αˆ=0.462,采用最大似然估计得到两参数的估计值分别为=1,αˆMLE=0.073 4；而运用本文所提出的贝叶斯估计为=0.972 7,αˆBayes=0.0718, 其中迭代次数N=1100, 舍弃迭代数N0=100, 起始点为α=1,σ=1, 先验选取无信息先验，可见后两种估计结果基本相当。

5 结束语

以上采用基于Gibbs抽样法获得了逐步增加II型截尾下Pareto分布的一种近似Bayes估计，可见Gibbs抽样是一种高效的马尔可夫链蒙特卡罗算法，其优点是避免了复杂的积分计算，缺点是获得的参数估计并不是精确的估计。

[1] 李凤, 师义民. 逐次定数截尾下Pareto分布参数的逆矩估计[J]. 统计与决策, 2010(24): 156-157.

[2] 李凤. 逐步增加Ⅱ型截尾下Pareto分布参数的区间估计[J]. 科学技术与工程, 2011,11(11): 2554-2557.

[3] WU S J, CHANG C T. Inference in the Pareto distribution based on progressive Type II censoring with random removals[J]. Applied Statistics, 2003, 30(2): 163-172.

[4] AMIN Z H. Bayesian inference for the Pareto lifetime model under progressive censoring with binomial removals[J]. Applied Statistics, 2008, 35(11): 1203-1217.

[5] SOLIMAN A A. Estimations for Pareto model using general progressive censored data and asymmetric loss[J]. Communications in Statistics—Theory and Methods, 2008, 37(9): 1353-1370.

[6] 武东, 汤银才. CE模型下指数分布步加试验的贝叶斯估计[J]. 数理统计与管理, 2008, 27(3): 423-427.

[7] BALAKRISHNAN N, AGGARWALA R. Progressive Censoring: Theory, Methods and Applications[M]. Boston: Birkhguser, 2000.

[8] NIGM A M, HAMDY H I. Bayesian prediction bounds for Pareto lifetime model[J]. Communications in Statistics—Theory and Methods,1987, 16(6): 1761-1772.

[9] 武东, 汤银才. 指数分布逐次定数截尾试验的多层贝叶斯估计[J]. 上海第二工业大学学报, 2011, 28(2): 113-116.

Bayesian Analysis of Pareto Distribution under Progressive Type-II Censoring

LI Qiong1, WU Dong2
( 1. Department of Electronic Information, Huishang Vocational Technical College, Hefei 230022, P. R. China; 2. School of Science, Anhui Agricultural University, Hefei 230036, P. R. China )

Bayesian analysis of Pareto distribution under ProgressiveType-II censoring was given. Bayesian estimation for parameters of the model was obtained using Markov chain Monte Carlo method. It is seen that Bayesian estimation is efficient through Monte Carlo simulation and an example.

progressive Type-II censoring; Pareto distribution; prior distribution; Bayesian estimation; Gibbs sampling

O213.2

A

1001-4543(2012)02-0135-04

2012-04-05;

2012-06-19

李琼（1983－），女，安徽泾县人，硕士，主要研究方向为贝叶斯统计、数据挖掘，电子邮箱qiongli1234@163.com。

国家自然科学基金资助项目（No. 10571057）；安徽省高等学校省级自然科学研究项目（No. KJ2011Z129）。