第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 Uk￾f 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] 将观测结束时 万方数据