第18卷第5期 系统仿真学报© Vol.18 No.5 2006年5月 Journal of System Simulation May,2006 基于MCMC稳态模拟的Weibul回归模型及其可靠性应用 林静，韩玉启'，朱慧明，陈杰1 (1.南京理工大学经济管理学院，江苏南京210094：2.湖南大学统计学系，湖南长沙410079 摘要：讨论了贝叶斯加速失效模型族中应用最为广泛的Weibull回归模型，提出针对寿命服从 Veibul分布的产品，运用基于Gibbs抽样的MCMC方法动态模拟出参数后验分布的马尔可夫链， 在失效率的先验分布为Gamma分布时，给出随机截尾条件下参数在Weibull回归模型中的贝叶斯 估计，提高了计算的精度。借助数据仿真说明了利用BUGS软件包进行建模分析的过程，证明了 该模型在可靠性应用中的直观性与有效性。 关键词：贝叶斯分析；可靠性；MCMC模拟；Gibbs抽样；Weibull分布 中图分类号：0212 文献标识码：A 文章编号：1004-731X(2006)05-1161-03 Weibull Regression Model Based on MCMC and Its Application in Reliability LIN Jing,HAN Yu-qi,ZHU Hui-ming,CHEN Jie (1.School of Economics Management,Nanjing University of Science Technology,Nanjing 210094,China; 2.Department of Statistics,Hunan University,Changsha 410079,China) Abstract:Weibull regression model was discussed,which is used widely in the family of Bayesian accelerated failure-time models.As for the productions whose life distributions belong to Weibull distribution,the MCMC method was brought forward based on Gibbs sampling to simulate dynamically the Markov Chain of the parameters'posterior distribution.From this,the parameters'Bavesian estimation of the Weibull regression model was given in the condition of the random truncated test and when the prior distribution of the failure rate belonged to the Gamma distribution.which improved the precision of the numeration.Also the data's simulation was utilized to show the process of setting the model by using the BUGS package It proves the objectivity and validity of the model. Key words:bayesian analysis;reliability;markov chain monte carlo simulation;gibbs sampling;weibull distribution 引言 为了避免高维数值积分，以往的许多研究都假设参数 为离散型随机变量，但这是一种很强的假设，它要求研究者 Weibull分布模型能够表现为各种不同的形状，较好的 必须对形状参数有较深入的认识：乔世君1在定数截尾试验 适用于各类数据，因而在可靠性分析中的应用十分广泛。实 中利用Q-Basic语言编写出贝叶斯估计的Gibbs抽样方案， 际中的许多场合需要我们确定产品寿命与某些主要伴随变 但仅仅得出有关参数的点估计，并没有得到更一般的结论： 量之间的关系，考察这种关系的方法之一是建立回归模型， 汤银才等利用数值积分技术得到Gibbs抽样的迭代结果， 在回归中使产品的寿命分布依赖于某些伴随变量（也称回归 但该方法在更高维条件下的计算仍然较为困难：韩明？提出 变量)。Weibull回归模型是加速失效模型族中非常重要的模 产品可靠度的E-Bayes估计，回避了数值积分问题，但其多 型，它将回归变量引入到寿命分布的描述，在应用中比单变 适用于无失效数据的情形。 量的Weibull模型更为有效。 随着计算机技术的发展和贝叶斯方法的改进，特别是马尔 近年来，为了适应小子样的情况以及充分利用已有的 可夫链蒙特卡罗Markov Chain Monte Carlo,,MCMC)方法I 验前信息，贝叶斯方法被广泛应用于可靠性参数的估计。特 以及BUGS(Bayesian inference Using Gibbs Sampling)软件 别是自Lindley和Smith四提出了多层先验分布的想法、韩 川的应用，原先异常复杂的数值计算问题游刃而解，参数后 明四提出先验分布的构造方法以来，贝叶斯方法在可靠性数 验分布的模拟也更为方便，现代贝叶斯理论及其应用日趋成 据的处理中取得了一些进展。然而，这些结果一般都要涉及 熟倒。本文在论述基于Gibbs抽样的MCMC方法的基础上， 数值积分技术。高维数值计算的困难，很大程度上制约了贝 借助数据仿真，给出了随机截尾条件下，利用BUGS软件 叶斯方法的应用。 包求解Weibull回归模型参数的贝叶斯估计的过程，提高了 收稿日期：2005-03-18 修回日期：2005-11-28 计算的精度。 基金项目：湖南省自然科学基金(05J0130) 作者简介：林静(1980-)，女，山东青岛人，博士生，研究方向为管 1 Weibull回归模型的贝叶斯分析 理科学与工程：韩玉启(1944)，男，江苏徐州人，教授、博导，研究 1.1随机截尾的Weibull回归模型 方向为管理方法论：朱慧明(1967-)，男，湖南湘潭人，博士后，研究 方向为贝叶斯统计应用：陈杰(1966-)，男，浙江泰化人，副教授，研 随机截尾试验例是一种简单且易于实现的研究过程：在 究方向为工业工程。 ·1161· 万方数据
第 18 卷第 5 期 系 统 仿 真 学 报© Vol. 18 No. 5 2006 年 5 月 Journal of System Simulation May 2006 • 1161 • 基于 MCMC 稳态模拟的 Weibull 回归模型及其可靠性应用 林 静 1 韩玉启 1 朱慧明 陈 杰 1 1.南京理工大学经济管理学院 江苏南京 210094 2.湖南大学统计学系 湖南长沙 410079 摘 要 讨论了贝叶斯加速失效模型族中应用最为广泛的 Weibull 回归模型 提出针对寿命服从 Weibull 分布的产品 运用基于 Gibbs 抽样的 MCMC 方法动态模拟出参数后验分布的马尔可夫链 在失效率的先验分布为 Gamma 分布时 给出随机截尾条件下参数在 Weibull 回归模型中的贝叶斯 估计 提高了计算的精度 借助数据仿真说明了利用 BUGS 软件包进行建模分析的过程 证明了 该模型在可靠性应用中的直观性与有效性 关键词 贝叶斯分析 可靠性 MCMC 模拟 Gibbs 抽样 Weibull 分布 中图分类号 O212 文献标识码 A 文章编号 1004-731X (2006) 05-1161-03 Weibull Regression Model Based on MCMC and Its Application in Reliability LIN Jing1 , HAN Yu-qi1 , ZHU Hui-ming2 , CHEN Jie1 (1.School of Economics & Management, Nanjing University of Science & Technology, Nanjing 210094, China; 2.Department of Statistics, Hunan University, Changsha 410079, China) Abstract: Weibull regression model was discussed, which is used widely in the family of Bayesian accelerated failure-time models. As for the productions whose life distributions belong to Weibull distribution, the MCMC method was brought forward based on Gibbs sampling to simulate dynamically the Markov Chain of the parameters’ posterior distribution. From this, the parameters’ Bayesian estimation of the Weibull regression model was given in the condition of the random truncated test and when the prior distribution of the failure rate belonged to the Gamma distribution, which improved the precision of the numeration. Also the data’s simulation was utilized to show the process of setting the model by using the BUGS package. It proves the objectivity and validity of the model. Key words: bayesian analysis; reliability; markov chain monte carlo simulation; gibbs sampling; weibull distribution 引 言1 Weibull 分布模型能够表现为各种不同的形状 较好的 适用于各类数据 因而在可靠性分析中的应用十分广泛 实 际中的许多场合需要我们确定产品寿命与某些主要伴随变 量之间的关系 考察这种关系的方法之一是建立回归模型 在回归中使产品的寿命分布依赖于某些伴随变量 也称回归 变量 Weibull 回归模型是加速失效模型族中非常重要的模 型 它将回归变量引入到寿命分布的描述 在应用中比单变 量的 Weibull 模型更为有效 近年来 为了适应小子样的情况以及充分利用已有的 验前信息 贝叶斯方法被广泛应用于可靠性参数的估计 特 别是自 Lindley 和 Smith [1] 提出了多层先验分布的想法 韩 明 [2] 提出先验分布的构造方法以来 贝叶斯方法在可靠性数 据的处理中取得了一些进展 然而 这些结果一般都要涉及 数值积分技术 高维数值计算的困难 很大程度上制约了贝 叶斯方法的应用 收稿日期 2005-03-18 修回日期 2005-11-28 基金项目 湖南省自然科学基金 05JJ0130 作者简介 林 静 1980- , 女, 山东青岛人, 博士生 研究方向为管 理科学与工程 韩玉启 1944- , 男, 江苏徐州人, 教授 博导 研究 方向为管理方法论 朱慧明 1967- , 男, 湖南湘潭人, 博士后 研究 方向为贝叶斯统计应用 陈 杰 1966- , 男, 浙江奉化人, 副教授 研 究方向为工业工程 为了避免高维数值积分 以往的许多研究都假设参数 为离散型随机变量 但这是一种很强的假设 它要求研究者 必须对形状参数有较深入的认识 乔世君 [3] 在定数截尾试验 中利用 Q-Basic 语言编写出贝叶斯估计的 Gibbs 抽样方案 但仅仅得出有关参数的点估计 并没有得到更一般的结论 汤银才 [ 4] 等利用数值积分技术得到 Gibbs 抽样的迭代结果 但该方法在更高维条件下的计算仍然较为困难 韩明 [5] 提出 产品可靠度的 E-Bayes 估计 回避了数值积分问题 但其多 适用于无失效数据的情形 随着计算机技术的发展和贝叶斯方法的改进 特别是马尔 可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)方法[ 6] 以及 BUGS (Bayesian inference Using Gibbs Sampling)软件 [7] 的应用 原先异常复杂的数值计算问题游刃而解 参数后 验分布的模拟也更为方便 现代贝叶斯理论及其应用日趋成 熟 [8] 本文在论述基于 Gibbs 抽样的 MCMC 方法的基础上 借助数据仿真 给出了随机截尾条件下 利用 BUGS 软件 包求解 Weibull 回归模型参数的贝叶斯估计的过程 提高了 计算的精度 1 Weibull 回归模型的贝叶斯分析 1.1 随机截尾的 Weibull 回归模型 随机截尾试验[9]是一种简单且易于实现的研究过程 在 万方数据
第18卷第5期 Vol.18 No.5 2006年5月 系统仿真学报 May,2006 n个样本中，设每个样本具有寿命T和截尾时间L,且T和 2马尔可夫链蒙特卡罗模拟 L是独立的连续型随机变量。假设试验中观测到产品寿命 y=minT,L},i=l,…,n,上式表明若T,>L,则第i个产 在贝叶斯分析中四，设k维随机向量U=(U,…,U) 品在L,时刻被截尾。规定v为指示变量，若T,≤L,v,=1, 具有联合分布π(U,,U),其中，U为模型参数或缺失 若T,>L,心，=0。则似然函数为： 的观测值，π)为后验分布。则对于我们感兴趣的函数(U) --Fn (1) 的数学期望为： EbU小-Jhu)z(u)du (6) 设某产品寿命y=(y,…y,)r服从二参数Weibull分布，记为 jπ(u)du W(a,y)。其概率密度为：fy/a,y)=a-exp(-y), 由于该积分往往形式复杂难于计算，此时采用蒙特卡罗积分 分布函数为：FUy/a,y)=1-exp(“)。为计算方便，令 进行近似，即： 元=logy),得到fUy/a,)=g-exp(1-exp(2)y),以 U小三U) (7) 及：Fy/a,)=1-exp(-exp(2)y)。则在随机截尾条件下， 当U,…,U相互独立时，由大数定律可知，样本容量n越大， 关于W(a,元)的似然函数如下： 其近似程度越高。但在很多复杂模型中，并不能简单地对 La.1/n.yv)If(y.la.2S(y.la.iy U,…,U,做出相互独立的假设，这就需要使用MCMC稳态 xa∑"exp{∑y2+ 模拟方法。MCMC模拟本质上是使用马尔可夫链的蒙特卡 罗积分，基本思想是：建立马尔可夫链对未知变量进行抽样 (.(e-Dlogtx,)-ex( (2) 模拟，当链达到稳态分布时即得所求的后验分布。基于贝叶 引入参数入的伴随变量即可建立加速失效模型族中十 斯推断原理的MCMC方法主要用于产生后验分布的样本， 分重要的模型Weibull回归模型：令入，=exp(xB),其中x, 计算边缘分布以及后验分布的矩。不同的抽样方法导致了不 为P×1维伴随变量，代表了影响样本寿命分布的主要因素； 同的MCMC方法，Gibbs抽样是其中最简单也是应用最广 B为p×1维回归系数向量，刻画了影响程度的大小，则 泛的一种。 Weibull回归模型的似然函数如下： Gibbs抽样过程属于马尔可夫更新机制的范畴。在上述 L(a,B/n,y,X,D)aexpxB+ v.(@-1)log(y.)-exp(xB)y"](3) 假设条件下，令U,代表某种随机变量或同组的几个随机变 量；第j组变量的边缘分布为(U)。给定任意初始向量 1.2回归参数的贝叶斯估计 Uo=(U,…,U),我们由fUIU,…,U)中抽取样本 当a己知时，exp(2)的共轭先验密度服从Gamma分 U;由fU,U",U,…,U)中抽取样本U四；由 布：当a与入均未知时，尚未有合适的联合先验密度， fU,1U,,U,U,…,U)抽取样本U";并由 此时，令a服从Gamma先验分布S(a。,ko),2服从正态先 fUIU",U,,U)抽取样本U":由上即完成了由 验分布N(4o,o),得到a与入的联合后验分布如下： U到U=(U",…,U)的转移。经过1次迭代，可以得到 π(a,元/m,D)xa,元/m,,u)x(a/(a/4,o) Uo=(U”,,U”),并最终得到U,U2,U,…。易证：由 =La.A/n.y.)( 'Ta)-e-ax 不同的U出发，当1→o,在遍历条件下，可以认为各时 1 exp--(2-BF) 刻U的边际分布为平稳分布，此时它收敛，并可以被看作 √2xo。 2 是样本的仿真观测点。而在收敛出现前的m次迭代中，各 xu-ep∑+2a-llog0W 状态的边际分布还不能认为是π(U),因此在估计Eh(U)】时 12-6P ep2）-xr-2 应将前m个迭代值去掉，即： (4) (8) 引入参数入的伴随变量，则Weibull回归模型a与B的联合 AUI立AU) n-m 后验分布如下： 3仿真分析 (a,B1m,y,X,)xa-xexp心，B+ 3.1随机截尾的Weibull回归模型 v,(a-1)log(y,)-exp(x;B)y)- 假定某种产品的寿命服从二参数Weibull分布W(a,Y), Ka-B-42(B-4) (5) 为了考察该产品在环境1与环境2中寿命的差异（即伴随变 量一环境对产品寿命的影响，且设仅有两种环境可能)，取 由上式可以看出，后验π(@，B/n,y,X,D)形式复杂，显然难 样本容量N=255,分别将其暴露于环境1和环境2中。在观 以给出其积分得精确形式。因而，本文对后验分布的计算借 助MCMC模拟方法。 察周期内记录个体失效的时间y,记为回：将观测结束时 ·1162· 万方数据
第 18 卷第 5 期 Vol. 18 No. 5 2006 年 5 月 系 统 仿 真 学 报 May, 2006 • 1162 • n 个样本中 设每个样本具有寿命T 和截尾时间 L 且T 和 L 是独立的连续型随机变量 假设试验中观测到产品寿命 min{ , }, i Ti Li y = i =1,L,n 上式表明若 , Ti > Li 则第 i 个产 品在 Li 时刻被截尾 规定u 为指示变量 若 , Ti £ Li ui = 1 若 , Ti > Li ui = 0 则似然函数为 i i f y F y n n L y n i n i i u u u - = = Õ - - å = 1 1 1 ( ) [1 ( )] ( )! ! ( ) 1 设某产品寿命 T 1 ( , ) n y = y y L 服从二参数 Weibull 分布 记为 W (a,g ) 其概率密度为 ( / , ) exp( ) a 1 a f y a g =agy -gy - 分布函数为 ( / , ) 1 exp( ) a F y a g = - gy 为计算方便 令 l = log(g ) 得到 ( / , ) exp( exp( ) ) a 1 a f y a l =ay l - l y - 以 及 ( / , ) 1 exp( exp( ) ) a F y a l = - - l y 则在随机截尾条件下 关于W (a,l) 的似然函数如下 1 { (1 ) 1 1 1 ( , /,, ) ( / , ) ( / , ) exp ( ( 1)log( ) exp( ) i i n i i n i i i n i i n i i i i L n y f y S y y y u u u a a l u a l a l a u l u a l = - = = = µ µ + å ü - - ý þ Õ å å 2 引入参数 l 的伴随变量即可建立加速失效模型族中十 分重要的模型 Weibull 回归模型 令 l exp( b) T i i = x 其中 i x 为 p ´1维伴随变量 代表了影响样本寿命分布的主要因素 b 为 p ´1维回归系数向量 刻画了影响程度的大小 则 Weibull 回归模型的似然函数如下 { ( 1)log( ) exp( ) ]} ( , / , , , ) exp [1 1 a u u a b a b u a u b i T i i i T i n i i y x y L n y X x n i i - - µ å = + å = 3 1.2 回归参数的贝叶斯估计 当a 已知时 exp(l) 的共轭先验密度服从 Gamma 分 布 当a 与 l 均未知时 尚未有合适的联合先验密度 [10] 此时 令a 服从 Gamma 先验分布 ( , ) z a0 k0 l 服从正态先 验分布 ( , ) 2 N m0 s 0 得到a 与 l 的联合后验分布如下 0 0 0 1 2 0 0 0 0 0 1 0 0 2 0 2 0 0 1 1 1 2 0 0 2 0 ( , /,, ) ( , /,,) ( / , ) ( / , ) ( , /,, ) ( exp( )) ( ) 1 1 ( exp( ( ) ) 2 2 exp ( ( 1)log( ) 1 exp( ) ) ( ) 2 n i i n n i i i i i i n y L n y L n y y y a a a u a p a l u a l u p aak p lms k a l u a k a a l m ps s a u l u a l k a l m s = - + - = = µ = × - ´ G - - å ì µ í + - - î ü - - - ý þ å å 4 引入参数 l 的伴随变量 则 Weibull 回归模型a 与 b 的联合 后验分布如下 0 1 1 1 1 0 0 0 0 ( , / , , , ) exp ( ( 1)log( ) exp( ) ) 1 ( ) ( ) 2 n i i n T i i i T i i i i T n y X x y x y a u a p a b u a u b u a b k a b m b m = + - = - ì µ å ´ + í î - - - ü - -S- ý þ å (5) 由上式可以看出 后验p (a,b / n, y, X,u) 形式复杂 显然难 以给出其积分得精确形式 因而 本文对后验分布的计算借 助 MCMC 模拟方法 2 马尔可夫链蒙特卡罗模拟 在贝叶斯分析中 [11] 设 k 维随机向量 ( , , ) U = U1 L Uk 具有联合分布 ( , , ) p U1 L Uk 其中 Uk 为模型参数或缺失 的观测值 p (×) 为后验分布 则对于我们感兴趣的函数 h(U ) 的数学期望为 [ ] ò ò = u du h u u du E h U ( ) ( ) ( ) ( ) p p 6 由于该积分往往形式复杂难于计算 此时采用蒙特卡罗积分 进行近似 即 [ ]» å = n t t h U n E h U 1 ( ) ( ) 1 ( ) 7 当U Uk , , 1 L 相互独立时 由大数定律可知 样本容量 n 越大 其近似程度越高 但在很多复杂模型中 并不能简单地对 U Uk , , 1 L 做出相互独立的假设 这就需要使用 MCMC 稳态 模拟方法 MCMC 模拟本质上是使用马尔可夫链的蒙特卡 罗积分 基本思想是 建立马尔可夫链对未知变量进行抽样 模拟 当链达到稳态分布时即得所求的后验分布 基于贝叶 斯推断原理的 MCMC 方法主要用于产生后验分布的样本 计算边缘分布以及后验分布的矩 不同的抽样方法导致了不 同的 MCMC 方法 Gibbs 抽样是其中最简单也是应用最广 泛的一种 Gibbs 抽样过程属于马尔可夫更新机制的范畴 在上述 假设条件下 令Ui 代表某种随机变量或同组的几个随机变 量 第 j 组变量的边缘分布为 ( ) U j f 给定任意初始向量 ( , , ) (0) (0) 1 (0) U = U L Uk 我们由 ( | , , ) (0) (0) U1 U2 Uk f L 中抽取样本 (1) U1 由 ( | , , , ) (0) (0) 3 (1) U2 U1 U Uk f L 中抽取样本 (1) U2 由 ( | , , (1) f U j U1 L , , , ) (0) (0) 1 (1) U j-1 U j+ L Uk 抽取样本 (1) U j 并 由 ( | , , , ) (1) 1 (1) 2 (1) Uk U1 U Ukf L 抽取样本 (1) Uk 由上即完成了由 (0) U 到 ( , , ) (1) (1) 1 (1) U = U L UK 的转移 经过 t 次迭代 可以得到 ( , , ) ( ) ( ) 1 ( ) t K t t U = U L U 并最终得到 U (1) ,U (2) ,U (3) ,L 易证 由 不同的 (0) U 出发 当 t ® ¥ 在遍历条件下 可以认为各时 刻 (t) U 的边际分布为平稳分布 此时它收敛 并可以被看作 是样本的仿真观测点 而在收敛出现前的 m 次迭代中 各 状态的边际分布还不能认为是p (U) 因此在估计 E[h(U)]时 应将前 m 个迭代值去掉 即 [ ] å - » = + n t m t h U n m E h U 1 ( ) ( ) 1 ( ) 8 3 仿真分析 3.1 随机截尾的 Weibull 回归模型 假定某种产品的寿命服从二参数 Weibull 分布W (a,g ) 为了考察该产品在环境 1 与环境 2 中寿命的差异 即伴随变 量 环境对产品寿命的影响 且设仅有两种环境可能 取 样本容量 N=255 分别将其暴露于环境 1 和环境 2 中 在观 察周期内记录个体失效的时间 i y 记为 t[i] 将观测结束时 万方数据
第18卷第5期 Vol.18 No.5 2006年5月 林静，等：基于MCMC稳态模拟的Weibull回归模型及其可靠性应用 May.2006 发生失效的时间视为有效数据，尚未失效的产品寿命是被随 信区间为(0.6911,0.7901)，B。的均值为-1.358,95%置 机截尾的，记为NA:只知道其寿命值不低于对其观测的时 信区间为(-1.621，-1.087)，B,的均值为-0.2635,95%置 间L,记为tcen[:伴随变量环境x的指示变量用trt[来 信区间为(-0.5823,0.04934)等。 表示，即当tt=0时表示产品样本暴露于环境1中，当tt=1 表2 WinBUGS中3000次抽样迭代的参数后验估计统计量 时表示暴露于环境2中（见表1）。因而，根据公式(2)得到 参数 均值 标准差 2.5% 中位数 97.5% B=(B。,B)',其中B。为截距，B为环境变量x,的系数： alpha 0.7925 0.05413 0.6911 0.7901 0.9063 设a与B得先验分布分别为：B~N,0,10), beta0 -1.358 0.1373 -1.621 -1.357 -1.087 a~S(L,0.00),即得一元Weibull回归模型。 betal -0.2635 0.1634 -0.5823 -0.2606 0.04934 表】随机截尾试验中产品寿命的仿真数据 alpha 样本i失效时间t0h裁尾时间tcen0h伴随指示变量t0 8.0 1 1.5 0.0 1 6.0 2 1.4 0.0 1 0 S55 3..... NA 7.3 ....… 0.8 D.7 0.8 beta0 254 0.9 0.0 0 4.0 255 NA 9.0 3.0 20 3.2基于BUGS的建模分析 2.0 -15 -10 BUGS是英国剑桥公共卫生研究所推出的利用MCMC方法 beta1 进行贝叶斯推断的专用软件包，使用BUGS可以很方便地 3.0 对许多常用的模型和分布进行Gibbs抽样，编程者只要设置 2.0 1.D 好变量的先验分布并对所研究的模型进行一般性的描述，就 00 能很容易实现对模型的贝叶斯分析。在BUGS中可以使用 10 -05 00 有向图模型方式(Directed Acyclic Graphical model,DAG)对 图2 WinBUGS中30O0次抽样迭代后参数的后验核密度估计 模型进行直观的描述（图1），也可以直接编写模型程序。 Gibbs抽样收敛后，可以得到参数后验分布的均值、标准差、 4结论 95%置信区间和中位数等信息（表2），并给出后验分布的核 Weibull回归模型在加速失效模型族中的应用十分广 密度估计图（图2）、参数的Gibbs抽样动态图、抽样值自相 泛，贝叶斯方法的应用，提高了该模型的有效性。以往的研 关图及均值和置信区间的变化图等，使抽样结果更直观、可 究由于不能很好的解决高维数值积分问题，从而制约了该理 靠：为了检验模型的收敛性，在BUGS中可以对有关参数 论在实际中的应用。本文利用基于Gibbs抽样的MCMC模 进行多条链的迭代分析，即同时输入多组初始值，经过 拟方法与BUGS软件的应用解决了该模型中高维数值计算 MCMC模拟形成多条迭代链，假如参数模型收敛（即符合 的不便，提高了计算的精度，有利于该模型在可靠性分析理 MCMC方法稳定条件时)则迭代图形的结果趋于重合。篇幅 论中的推广。本文仿真实例中的一元Weibull回归模型可以 所限，文中仅列出其中一部分。这里基于表1仿真数据，建 推广到多元的情况中去，不难看出，该模型在可靠性分析领 立如上的一元Weibull回归模型，截取前1000次迭代结果， 域中的应用前景比单纯Weibull模型更为广泛。需要指出， 从第1001次开始进行3000次迭代分析。 在仿真分析中，我们给定的先验分布的参数形式也具有普遍 由BUGS运行结果可以得到a的均值为0.7925,95%置 性，实际中应具体情况具体分析。 参考文献： [1]Lindley D V,Smith A F M.Bayes estimation for the linear model [J]. Journal of the Royal Statistical Society Series B(1369-7412),1972 34:1-41. [2]韩明.多层先验分布的构造及其应用[U运筹与管理，1997,6(3) 31-40. [3)乔世君，张世英.用Gibbs抽样算法计算定数截尾时Weibull分布 的贝叶斯估计[).数理统计与管理，2000,192)：35-40. for(IN1:N) [4)汤银才，费鹤良.基于Gibbs抽样的Weibull分布序进应力加速寿 命试验的Bayes分析.数理统计与应用概率，1998,13(1)：81-88 图1 WinBUGS中的贝叶斯有向无环图 ·1163· (下转第1185页) 万方数据
第18卷第5期 Vol. 18 No. 5 2006 年 5 月 林 静, 等 基于 MCMC 稳态模拟的 Weibull 回归模型及其可靠性应用 May, 2006 • 1163 • 发生失效的时间视为有效数据 尚未失效的产品寿命是被随 机截尾的 记为 NA 只知道其寿命值不低于对其观测的时 间 Li 记为 t.cen[i] 伴随变量环境 x 的指示变量用 trt[i]来 表示 即当 trt=0 时表示产品样本暴露于环境 1 中 当 trt=1 时表示暴露于环境 2 中(见表 1) 因而 根据公式(2)得到 T ( , ) = bbb 10 其中 b0 为截距 b1 为环境变量 i x 的系数 设 a 与 b 得先验分布分别为 ~ (0,10 ) 4 2 b N I a ~ z (1,0.001) 即得一元 Weibull 回归模型 表 1 随机截尾试验中产品寿命的仿真数据 样本i 失效时间t[i]/h 截尾时间t.cen[i]/h 伴随指示变量trt[i] 1 2 3 M M 254 255 1.5 1.4 NA M M 0.9 NA 0.0 0.0 7.3 M M 0.0 9.0 1 1 0 M M 0 1 3.2 基于 BUGS 的建模分析 BUGS 是英国剑桥公共卫生研究所推出的利用 MCMC 方法 进行贝叶斯推断的专用软件包 使用 BUGS 可以很方便地 对许多常用的模型和分布进行 Gibbs 抽样 编程者只要设置 好变量的先验分布并对所研究的模型进行一般性的描述 就 能很容易实现对模型的贝叶斯分析 在 BUGS 中可以使用 有向图模型方式(Directed Acyclic Graphical model, DAG )对 模型进行直观的描述 图 1 也可以直接编写模型程序 Gibbs 抽样收敛后 可以得到参数后验分布的均值 标准差 95%置信区间和中位数等信息 表 2 并给出后验分布的核 密度估计图 图 2 参数的 Gibbs 抽样动态图 抽样值自相 关图及均值和置信区间的变化图等 使抽样结果更直观 可 靠 为了检验模型的收敛性 在 BUGS 中可以对有关参数 进行多条链的迭代分析 即同时输入多组初始值 经过 MCMC 模拟形成多条迭代链 假如参数模型收敛(即符合 MCMC 方法稳定条件时)则迭代图形的结果趋于重合 篇幅 所限 文中仅列出其中一部分 这里基于表 1 仿真数据 建 立如上的一元 Weibull 回归模型 截取前 1000 次迭代结果 从第 1001 次开始进行 3000 次迭代分析 由 BUGS 运行结果可以得到a 的均值为 0.7925 95%置 图 1 WinBUGS 中的贝叶斯有向无环图 信区间为 0.6911 0.7901 b0 的均值为-1.358 95%置 信区间为 -1.621 -1.087 b1 的均值为-0.2635 95%置 信区间为 -0.5823 0.04934 等 表 2 WinBUGS 中 3000 次抽样迭代的参数后验估计统计量 参数 均值 标准差 2.5% 中位数 97.5% alpha beta0 beta1 0.7925 -1.358 -0.2635 0.05413 0.1373 0.1634 0.6911 -1.621 -0.5823 0.7901 -1.357 -0.2606 0.9063 -1.087 0.04934 图 2 WinBUGS 中 3000 次抽样迭代后参数的后验核密度估计 4 结论 Weibull 回归模型在加速失效模型族中的应用十分广 泛 贝叶斯方法的应用 提高了该模型的有效性 以往的研 究由于不能很好的解决高维数值积分问题 从而制约了该理 论在实际中的应用 本文利用基于 Gibbs 抽样的 MCMC 模 拟方法与 BUGS 软件的应用解决了该模型中高维数值计算 的不便 提高了计算的精度 有利于该模型在可靠性分析理 论中的推广 本文仿真实例中的一元 Weibull 回归模型可以 推广到多元的情况中去 不难看出 该模型在可靠性分析领 域中的应用前景比单纯 Weibull 模型更为广泛 需要指出 在仿真分析中 我们给定的先验分布的参数形式也具有普遍 性 实际中应具体情况具体分析 参考文献: [1] Lindley D V, Smith A F M. Bayes estimation for the linear model [J]. Journal of the Royal Statistical Society Series B(1369-7412), 1972, 34: 1-41. [2] 韩明.多层先验分布的构造及其应用[J].运筹与管理, 1997, 6(3): 31-40. [3] 乔世君, 张世英. 用 Gibbs 抽样算法计算定数截尾时 Weibull 分布 的贝叶斯估计[J]. 数理统计与管理, 2000, 19(2): 35-40. [4] 汤银才 费鹤良. 基于 Gibbs 抽样的 Weibull 分布序进应力加速寿 命试验的 Bayes 分析[J]. 数理统计与应用概率, 1998, 13(1):81-88. 下转第 1185 页 万方数据
第18卷第5期 Vol.18 No.5 2006年5月 李振宇，等：一种结构化P2P系统的拓扑匹配算法 May,2006 较多的RTT探测，所以在这种拓扑中，算法的开销较大。 15 4结论 -M40 5 1 -*2700 本文提出了一种结构化P2P系统中逻辑覆盖网络和物 理网络之间拓扑匹配算法：事件驱动算法。仿真实验表明事 05 件驱动算法可以在开销很小的情况下使系统的逻辑链路延 迟和物理链路延迟的比值降低65%左右。与其他现有算法 2040608000120 相比，事件驱动算法优势在于：1)吸收了目前存在的两种 Time in Minute 思路的优点，克服了每种思路的不足，从而大大提高了算法 图8物理网络规模对算法开销的影响 的有效性和准确性：2)算法带来的额外开销是非常小的。 E。物理网络拓扑的类型对算法的影响 进一步的研究将主要集中在非结构化P2P系统的拓扑 图9和图10中系统参数为d=4,lee=75,join=75。 匹配，在此基础上将对资源的快速定位算法做深入的研究。 25 参考文献： [1]Napster[EB/OL].http://www.napster.com/. [2]Gnutella[EB/OL].http://www.gnutella.com/. [3]S Jiang,L Guo,X Zhang.Ligh ood:an ef cient ooding scheme for le search in unstructured peer-to-peer systems[Cl/In Proceedings of -4-t3k-smal 13k-lorge ICPP 2003,Kaohsiung,Taiwan,October 2003:149-160. [4]S Ratnasamy,P Francis,M Handley,R Karp.A scalable content-addressable network[CV/In Proceedings of SIGCOMM 2001. 20 9nmne时 100120 San Diego,CA,USA,August 2001:161-172. 图9物理网络类型对stretch的影响 [5]I Stoica,R Morris,D Karger,M Kaashoek,H Balakrishnan.Chord:A scalable peer-to-peer lookup service for internet applications[Cl/.In Proceedings of SIGCOMM 2001.San Deigo,CA,USA,August 2001: 149-160. [6]M.Castro,P.Druschel,Y.Hu,and A.Rowstron.Exploiting network proximity in distributed hash tables[C/In Proceedings of FuDiCo 2002,Bertinoro,Italy,June 2002 -。t3k-5nl [7]B Zhao,J Kubiatowicz,A Joseph.Tapestry:An infrastructure for +t3米93 fault-tolerant wide-area location and routing[R.Technical report,UC Berkeley,April 2001. 20 40 100 120 [8]S Ratnasamy,M Handley,R Karp,S Shenker.Topologically-aware 图10物理网络类型对算法开销的影响 overlay construction and server selection[Cl/In Proceedings of INFOCOM 2002.New York,NY,USA,June 2002. 从图9和图I0可以看出，就系统的stretch而言， [9]Z Xu,C Tang,Z Zhang.Building topology-aware overlays using t3k-small优于t3k-large:而就算法开销而言，t3k-small比 global soft-state[Cl/In Proceedings of ICDCS 2003.Providence,RI. t3k-large略大。前面己经提到在t3k-small中，每个stub域 USA,May2003:500-508. [10]Shansi Ren,Lei Guo,Song Jiang.Xiaodong Zhang.SAT-Match:a 内的节点个数比较多，而这些节点间的延迟相差很小，所以 self-adaptive topology matching method to achieve low lookup latency 一个sub域中的节点落在同一个界标区域的概率比较大。 in structured P2P overlay networks[Cl/.IPDPS'04.April 26-30,2004. 这就会使一个节点在有限的探测内可能找到距离该节点最 [11]E W Zegura,K L Calvert,S Bhattacharjee.How to model an 近的节点，从而有利于系统stretch的降低。在t3k-small中， internetwork[Cl/In Proceedings of INFOCOM 1996.San Francisco. CA,USA,March 1996,volume 2,594-602 一个stub域中的节点以很高的概率落入同一个界标区域中， 而且每个stub域中的节点个数比较多，这就导致需要做比 一命一一一一一一一一一一一一一一一 一一一一一一一一一一一一一一一 (上接第1163页) 济管理，2004，(6)：50-54. [)】韩明.可靠度的一种新估计方法.兵工学报.2004,25(1上60-64. [9例Lawless J F.寿命数据中的统计模型与方法M.茆诗松等译.北 [6]Gilks W R,Richardson S,Spiegelhalter D J.Markov Chain Monte 京：中国统计出版社，1998. Carlo in Practice [M].London:Chapman Hall,1996. [10]Perter Congdon.Applied Bayesian Modeling [M].England:John [7]David S.Andrew T,Nicky B.BUGS.5:Bayesian inference Using Wiley Sons Ltd,2003. Gibbs Sampling Manual (version ii)[M].Cambridge,U.K.:MRC [11]Ibrahim J G,Chen M H,Sinha D.Bayesian Survival Analysis [M] Biostatistics Unit,1996 刘乐平，袁卫.现代贝叶斯分析与现代统计推断.经济理论与经 New York Berlin Heidelberg,2001 ·1185· 万方数据
第18卷第5期 Vol. 18 No. 5 2006 年 5 月 李振宇, 等 一种结构化 P2P 系统的拓扑匹配算法 May, 2006 • 1185 • 图 8 物理网络规模对算法开销的影响 E. 物理网络拓扑的类型对算法的影响 图 9 和图 10 中系统参数为 d=4, leave=75 join=75 图 9 物理网络类型对 stretch 的影响 图 10 物理网络类型对算法开销的影响 从图 9 和图 10 可以看出 就系统的 stretch 而言 t3k-small 优于 t3k-large 而就算法开销而言 t3k-small 比 t3k-large 略大 前面已经提到在 t3k-small 中 每个 stub 域 内的节点个数比较多 而这些节点间的延迟相差很小 所以 一个 stub 域中的节点落在同一个界标区域的概率比较大 这就会使一个节点在有限的探测内可能找到距离该节点最 近的节点 从而有利于系统 stretch 的降低 在 t3k-small 中 一个 stub 域中的节点以很高的概率落入同一个界标区域中 而且每个 stub 域中的节点个数比较多 这就导致需要做比 较多的 RTT 探测 所以在这种拓扑中 算法的开销较大 4 结论 本文提出了一种结构化 P2P 系统中逻辑覆盖网络和物 理网络之间拓扑匹配算法 事件驱动算法 仿真实验表明事 件驱动算法可以在开销很小的情况下使系统的逻辑链路延 迟和物理链路延迟的比值降低 65 左右 与其他现有算法 相比 事件驱动算法优势在于 1 吸收了目前存在的两种 思路的优点 克服了每种思路的不足 从而大大提高了算法 的有效性和准确性 2 算法带来的额外开销是非常小的 进一步的研究将主要集中在非结构化 P2P 系统的拓扑 匹配 在此基础上将对资源的快速定位算法做深入的研究 参考文献 [1] Napster[EB/OL]. http://www.napster.com/. [2] Gnutella[EB/OL]. http://www.gnutella.com/. [3] S Jiang, L Guo, X Zhang. Ligh ood: an ef cient ooding scheme for le search in unstructured peer-to-peer systems[C]// In Proceedings of ICPP 2003, Kaohsiung, Taiwan, October 2003: 149–160. [4] S Ratnasamy, P Francis, M Handley, R Karp. A scalable content-addressable network[C]// In Proceedings of SIGCOMM 2001. San Diego, CA, USA, August 2001: 161–172. [5] I Stoica, R Morris, D Karger, M Kaashoek, H Balakrishnan. Chord: A scalable peer-to-peer lookup service for internet applications[C]//. In Proceedings of SIGCOMM 2001. San Deigo, CA, USA, August 2001: 149–160. [6] M. Castro, P. Druschel, Y. Hu, and A. Rowstron. Exploiting network proximity in distributed hash tables[C]// In Proceedings of FuDiCo 2002, Bertinoro, Italy, June 2002 [7] B Zhao, J Kubiatowicz, A Joseph. Tapestry: An infrastructure for fault-tolerant wide-area location and routing[R]. Technical report, UC Berkeley, April 2001. [8] S Ratnasamy, M Handley, R Karp, S Shenker. Topologically-aware overlay construction and server selection[C]// In Proceedings of INFOCOM 2002. New York, NY, USA, June 2002. [9] Z Xu, C Tang, Z Zhang. Building topology-aware overlays using global soft-state[C]// In Proceedings of ICDCS 2003. Providence, RI, USA, May 2003: 500–508. [10] Shansi Ren, Lei Guo, Song Jiang, Xiaodong Zhang. SAT-Match: a self-adaptive topology matching method to achieve low lookup latency in structured P2P overlay networks[C]//. IPDPS'04. April 26-30, 2004. [11] E W Zegura, K L Calvert, S Bhattacharjee. How to model an internetwork[C]// In Proceedings of INFOCOM 1996. San Francisco, CA, USA, March 1996, volume 2, 594–602. 上接第 1163 页 [5] 韩明. 可靠度的一种新估计方法[J]. 兵工学报. 2004, 25(1): 60-64. [6] Gilks W R, Richardson S, Spiegelhalter D J. Markov Chain Monte Carlo in Practice [M]. London: Chapman & Hall, 1996. [7] David S, Andrew T, Nicky B. BUGS 0.5: Bayesian inference Using Gibbs Sampling Manual (version ii)[M]. Cambridge, U.K.: MRC Biostatistics Unit, 1996. [8] 刘乐平,袁卫. 现代贝叶斯分析与现代统计推断[J]. 经济理论与经 济管理, 2004, (6): 50-54. [9] Lawless J F. 寿命数据中的统计模型与方法[M]. 茆诗松等译. 北 京: 中国统计出版社, 1998. [10] Perter Congdon. Applied Bayesian Modeling [M]. England: John Wiley & Sons Ltd, 2003. [11] Ibrahim J G, Chen M H, Sinha D. Bayesian Survival Analysis [M]. New York Berlin Heidelberg, 2001. 万方数据
