第28卷第3期 兵工 学报 Vol.28 No.3 2007年3月 ACTA ARMAMENTARII Mar.2007 基于MCMC稳态模拟方法的弹药贮存可靠性评估模型 林静1，韩玉启1，朱慧明2，王晔3 (1.南京理工大学，江苏南京210094：2.湖南大学，湖南长沙410079： 3.海军工程大学，湖北武汉430033) 摘要：基于贝叶斯生存分析理论的日渐成熟，以及传统弹药贮存可靠性评估方法的不足，提出 在弹药贮存可靠性的随机截尾试验中，引入贝叶斯生存分析的相关理论，构建出可靠度的回归模 型。以常用的指数分布与Weibull分布为例，运用基于Gibbs抽样的MCMC方法动态模拟出参数 后验分布的马尔可夫链，给出随机截尾条件下参数的贝叶斯估计。利用BUGS软件包进行建模仿 真分析的结果，证明该模型在弹药贮存可靠性评估中的直观性与有效性。 关键词：应用统计数学：贝叶斯分析：弹药：可靠性：MCMC方法：贮存 中图分类号：O212 文献标志码：A 文章编号：1000-1093(2007)03-0315-04 Reliability Models of Ammunition Storage Based on MCMC Simulation Method LIN Jing',HAN Yu-qi,ZHU Hui-ming2,WANG Ye3 (1.Nanjing University of Science Technology,Nanjing 210094,Jiangsu,China: 2.Hunan University,Changsha 410079,Hunan,China: 3.Naval University of Engineering,Wuhan 430033,Hubei,China) Abstract:With the more and more development of the Bayesian survival analysis theory,as well as the fault of the traditional methods for storage reliability evaluation of ammunition,the Bayesian analysis method was proposed to build reliability regression models in the condition of the random truncated test.Take exponential distribution and Weibull distribution for example,and the MCMC method based on Gibbs sampling was put forward to simulate dynamically the Markov chain of the parameters' posterior distribution.Also,the parameters'Bayesian estimations were given out.The BUGS package was used to realize the simulation.The results prove the model to be objective and valid in the storage reliability evaluation of ammunition. Key words:applied statistical mathematics:Bayesian analysis:ammunition:reliability:MCMC method:storage 为了随时监控弹药贮存的可靠性，文献[1]首次 环境、阵地环境等)下弹药贮存可靠性进行评估时， 提出运用能够有效处理各类删失数据的生存分析方 大都采用相对独立的评估模型，而鲜有对同种弹药 法，使评估结果更加接近产品实际：进而文献[2]更 在不同各种环境因子下寿命异同的集成评估3-4]。 为系统的论证了生存分析法在弹药贮存可靠性评估 随着计算机技术的发展和贝叶斯方法的改进（主要 中解决弹药寿命数据删失的应用。然而，近年来，为 是马尔可夫链蒙特卡罗(Markov Chain Monte Car- 了适应小子样的情况以及充分利用己有的验前信 lo,MCMC)方法以及BUGS(Bayesian inference Us- 息，贝叶斯方法与生存分析相结合产生了贝叶斯生 ing Gibbs Sampling)软件的应用，原先异常复杂的数 存分析理论，并已成为国内外统计研究的新热点：另 值计算问题游刃而解：特别是在需要建立回归模型 一方面，国内外学者在对不同环境因子（如各种自然 以考察系统寿命与某些伴随变量（也称回归变量）之 收稿日期：2006-04-18 基金项目：湖南省自然科学基金(05J0130)
基于 MCMC 稳态模拟方法的弹药贮存可靠性评估模型 林静1 ，韩玉启1 ，朱慧明2 ，王晔3 （1. 南京理工大学，江苏 南京210094 ；2. 湖南大学，湖南 长沙410079 ； 3. 海军工程大学，湖北 武汉430033 ） 摘要：基于贝叶斯生存分析理论的日渐成熟，以及传统弹药贮存可靠性评估方法的不足，提出 在弹药贮存可靠性的随机截尾试验中，引入贝叶斯生存分析的相关理论，构建出可靠度的回归模 型。以常用的指数分布与 Wei bull 分布为例，运用基于 Gi bbs 抽样的 MCMC 方法动态模拟出参数 后验分布的马尔可夫链，给出随机截尾条件下参数的贝叶斯估计。利用BUGS 软件包进行建模仿 真分析的结果，证明该模型在弹药贮存可靠性评估中的直观性与有效性。 关键词：应用统计数学；贝叶斯分析；弹药；可靠性；MCMC 方法；贮存 中图分类号：O212 文献标志码：A 文章编号：1000-1093（2007 ）03-0315-04 Reliabilit y Models of ammunition storage based on MCMC si mulation Met hod LI N Ji ng1 ，~AN Yu-Ci 1 ，Z~U ~ui- mi ng2 ，WANG Ye3 （1. Nanji ng University of Science & technology ，Nanji ng 210094 ，Jiangsu ，Chi na ； 2. ~unan University ，Changsha 410079 ，~unan ，Chi na ； 3. Naval University of engi neeri ng ，Wuhan 430033 ，~ubei ，Chi na ） Abstr act ：Wit h t he more and more develop ment of t he Bayesian survi val analysis t heory ，as well as t he f ault of t he traditional met hods f or storage reliability eval uation of a mmunition ，t he Bayesian analysis met hod was proposed to build reliability regression models i n t he condition of t he random truncated test . take exponential distri bution and Wei bull distri bution f or exa mple ，and t he MCMC met hod based on Gi bbs sa mpli ng was put f or ward to si mulate dyna mically t he Markov chai n of t he para meters ’ posterior distri bution . Also ，t he para meters ’Bayesian esti mations were gi ven out . the BUGSpackage was used to realize t he si mulation . the results prove t he model to be obj ecti ve and vali d i n t he storage reliability eval uation of a mmunition . Key Wor ds ：applied statistical mat hematics ；Bayesian analysis ；a mmunition ；reliability ；MCMC met hod ；storage 收稿日期：2006 -04 -18 基金项目：湖南省自然科学基金（05JJ0130 ） 为了随时监控弹药贮存的可靠性，文献［1 ］首次 提出运用能够有效处理各类删失数据的生存分析方 法，使评估结果更加接近产品实际；进而文献［2 ］更 为系统的论证了生存分析法在弹药贮存可靠性评估 中解决弹药寿命数据删失的应用。然而，近年来，为 了适应小子样的情况以及充分利用已有的验前信 息，贝叶斯方法与生存分析相结合产生了贝叶斯生 存分析理论，并已成为国内外统计研究的新热点；另 一方面，国内外学者在对不同环境因子（如各种自然 环境、阵地环境等）下弹药贮存可靠性进行评估时， 大都采用相对独立的评估模型，而鲜有对同种弹药 在不同各种环境因子下寿命异同的集成评估［3 -4 ］。 随着计算机技术的发展和贝叶斯方法的改进（主要 是马尔可夫链蒙特卡罗（Markov Chai n Monte Carlo ，MCMC）方法以及BUGS（Bayesian i nf erence Usi ng Gi bbs Sa mpli ng ）软件的应用，原先异常复杂的数 值计算问题游刃而解；特别是在需要建立回归模型 以考察系统寿命与某些伴随变量（也称回归变量）之 第28 卷第3 期 2 0 0 7 年3 月 兵 工 学 报 ACtA ARMAMeNtARII Vol .28 No .3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Mar . 2007
316 兵工学报 第28卷 间关系的场合，参数后验分布的模拟更为方便，贝叶 而在收敛出现前的m次迭代中，各状态的边际分布 斯生存分析理论及其应用日趋成熟。 还不能认为是π(U),因此在估计E[h(U)]时应将 在此背景下，本文提出利用贝叶斯生存分析的 前m个迭代值去掉，即： 有关理论，针对文献[1]中的模型从2个方面加以改 进：1)在模型中引入基于Gibbs抽样的MCMC稳 E[h(U)]≈1 >】h(Ux). (3) n-m.=m+1 态模拟方法，构建弹药贮存可靠性的贝叶斯生存分 析模型，以解决传统贝叶斯分析中高维数值积分的 2指数回归模型的贝叶斯分析 困难：2)构建弹药贮存可靠度的回归模型，有效地 2.1随机截尾的指数回归模型 反映出不同环境条件对弹药贮存可靠性的影响。 文献[1]中指出，弹药贮存可靠性试验的一种常 1基于Gibbs的MCMC模拟方法[5-6] 见情况是随机试验。在该类试验中，发现弹药失效 的时刻并非弹药实际失效的时刻，弹药的寿命数据 设k维随机向量U=(U1,…,U)具有联合分 具有截尾的特性，即：只知道弹药个体中一部分寿命 布π(U1,…,U5),其中，U(i=1,…,k)为模型参 值低于某个值，而剩余部分的寿命只知其超过某一 数或缺失的观测值，π(·)为其后验分布，函数h(U) 特定值。截尾方式的不同导致了不同的寿命试验及 的数学期望为： 研究方法的各异]。弹药贮存的随机截尾试验是 很简单且经常可实现的一种研究过程：将试验中发 E[h(U)]=h(u)x(u)du x(u)du, 现弹药失效的年份减1，规定为该枚弹药的寿命值， (1) 则在枚弹药样本中，每个样本具有寿命Y和截 由于该积分往往形式复杂难于计算，此时采用 尾时间L(如：在第5年试验失效的弹药规定其寿命 蒙特卡罗积分进行近似，即： Y=4,若该枚弹药未失效，则具有截尾寿命时间 ELh(U)]之 L=5),Y和L是独立的连续随机变量。设试验中 h(U) (2) 能且仅能观测到弹药样本的寿命t:=minY;,L,J, 当U1,…,U相互独立时，由大数定律可知， (i=1,…,n).若Y,>L;,则第i枚弹药在L,后失 样本容量越大，其近似程度越高。但在很多复杂 去观察（被截尾）。定义指示变量：当Y:≤L= 模型中，并不能简单地对U1,…,U做出相互独立 1,当Y;>L,=0.于是得到弹药贮存寿命在随机 的假设，这就需要使用MCMC稳态模拟方法。 截尾试验中的似然函数为： MCMC模拟本质上是使用马尔可夫链的蒙特卡罗 L(t)= n f(t)[1-F(t)]-". 积分，基本思想是：建立马尔可夫链对未知变量进行 抽样模拟，当链达到稳态分布时即得所求的后验分 i=1 布。基于贝叶斯推断原理的MCMC方法主要用于 (4) 产生后验分布的样本，计算边缘分布以及后验分布 对于文献[1]中所指出的、贮存寿命服从指数分 的矩。不同的抽样方法导致了不同的MCMC方法， 布的弹药，设其贮存寿命t=（t1,…,tn)T服从失效 Gibbs抽样是其中最简单也是应用最广泛的一种。 率为入的指数分布，记为e(λ)，则其失效时间密度 Gibbs抽样过程属于马尔可夫更新机制的范 f(t:|入)=λexp(-t入)，可靠度及其失效时间分布 畴。在上述假设条件下，令U,代表某种随机变量 函数为：F(t:|a)=1-exp(-t入)=1-R(t|入). 或同组的几个随机变量：第；组变量的边缘分布为 综上可得随机截尾条件下，关于ε（入）的似然函数如 f(U).给定任意初始向量U0)=(U0),, 下： U9),由f(U1lU,…,U9)中抽取样本U: L(ln,t,)II ftl)"R( 由f(U2UP,U,,U0)中抽取样本U”:由 i=1 f(UU,…,U,U91,…,U9)抽取样本 Ⅱ[，exp(-t)][exp(-ta)J1-.(5) U:并由f(UsU,U,,U1)抽取样本 文献[1]中将贮存可靠度表示为故有可靠度的 U:由上即完成了由Uo到U1)=(U,…, 条件贮存可靠度的乘积，认为固有可靠度应由弹药 U)的转移。经过之次迭代，可以得到U)= 出厂验收的试验数据得到。实际中，由于贮存环境 (U,…,U),并最终得到U1,U2),U3”,… 的复杂性使得该方法的局限性极强。由此引入失效 易证[5：由不同的U0出发，当之→∞，在遍历条件 率入的伴随变量，建立指数回归模型，以确定贮存 下，可以认为各时刻Ux)的边际分布为平稳分布， 寿命与某些环境变量（如：湿度、温度、辐射等）之间 此时它收敛，并可以被看作是样本的仿真观测点。 的关系：令λ；=exp(xB),式中：x:为伴随变量，代
间关系的场合，参数后验分布的模拟更为方便，贝叶 斯生存分析理论及其应用日趋成熟。 在此背景下，本文提出利用贝叶斯生存分析的 有关理论，针对文献［1 ］中的模型从2 个方面加以改 进：1 ）在模型中引入基于 Gi bbs 抽样的 MCMC 稳 态模拟方法，构建弹药贮存可靠性的贝叶斯生存分 析模型，以解决传统贝叶斯分析中高维数值积分的 困难；2 ）构建弹药贮存可靠度的回归模型，有效地 反映出不同环境条件对弹药贮存可靠性的影响。 1 基于Gibbs的 MCMC模拟方法［5 -6 ］ 设a 维随机向量!=（U1 ，…，Ua ）具有联合分 布!（U1 ，…，Ua ），其中，Ui （i = 1 ，…，a ）为模型参 数或缺失的观测值，!（·）为其后验分布，函数h（!） 的数学期望为： E［h（!）］=｝h（ ）!（ ）d （｝!（ ）d ）， （1 ） 由于该积分往往形式复杂难于计算，此时采用 蒙特卡罗积分进行近似，即： E［h（!）］＞1 7 】 7 Z =1 h（! （Z ））. （2 ） 当 U1 ，…，Ua 相互独立时，由大数定律可知， 样本容量7 越大，其近似程度越高。但在很多复杂 模型中，并不能简单地对 U1 ，…，Ua 做出相互独立 的假 设，这 就 需 要 使 用 MCMC 稳 态 模 拟 方 法。 MCMC 模拟本质上是使用马尔可夫链的蒙特卡罗 积分，基本思想是：建立马尔可夫链对未知变量进行 抽样模拟，当链达到稳态分布时即得所求的后验分 布。基于贝叶斯推断原理的 MCMC 方法主要用于 产生后验分布的样本，计算边缘分布以及后验分布 的矩。不同的抽样方法导致了不同的 MCMC 方法， Gi bbs 抽样是其中最简单也是应用最广泛的一种。 Gi bbs 抽样过程属于马尔可夫更新机制的范 畴。在上述假设条件下，令 Ui 代表某种随机变量 或同组的几个随机变量；第 组变量的边缘分布为 （U ）. 给 定 任 意 初 始 向 量 ! （0 ）=（U（0 ） 1 ，…， U（0 ） a ），由 （U1】U（0 ） 2 ，…，U（0 ） a ）中抽取样本 U（1 ） 1 ； 由 （U2】U（1 ） 1 ，U（0 ） 3 ，…，U（0 ） a ）中抽取样本 U（1 ） 2 ；由 （U 】U（1 ） 1 ，…，U（1 ） -1 ，U（0 ） + 1 ，…，U（0 ） a ）抽取样本 U（1 ） ；并由 （Ua】U（1 ） 1 ，U（1 ） 2 ，…，U（1 ） a -1 ）抽取样本 U（1 ） a ；由上即完成了由 ! （0 ）到 ! （1 ） =（U（1 ） 1 ，…， U（1 ） a ）的转移。经过 Z 次迭代，可以得到 ! （Z ） = （U（Z ） 1 ，…，U（Z ） a ），并最终得到 ! （1 ），! （2 ），! （3 ），…. 易证［5 ］：由不同的 ! （0 ）出发，当Z 一 ，在遍历条件 下，可以认为各时刻 ! （Z ）的边际分布为平稳分布， 此时它收敛，并可以被看作是样本的仿真观测点。 而在收敛出现前的 m 次迭代中，各状态的边际分布 还不能认为是!（!），因此在估计E［h（!）］时应将 前 m 个迭代值去掉，即： E［h（!）］＞ 1 7 - m 】 7 Z = m+1 h（! （Z ））. （3 ） 2 指数回归模型的贝叶斯分析 2.1 随机截尾的指数回归模型 文献［1 ］中指出，弹药贮存可靠性试验的一种常 见情况是随机试验。在该类试验中，发现弹药失效 的时刻并非弹药实际失效的时刻，弹药的寿命数据 具有截尾的特性，即：只知道弹药个体中一部分寿命 值低于某个值，而剩余部分的寿命只知其超过某一 特定值。截尾方式的不同导致了不同的寿命试验及 研究方法的各异［7 ］。弹药贮存的随机截尾试验是 很简单且经常可实现的一种研究过程：将试验中发 现弹药失效的年份减1 ，规定为该枚弹药的寿命值， 则在7 枚弹药样本中，每个样本具有寿命 Y 和截 尾时间L（如：在第5 年试验失效的弹药规定其寿命 Y = 4 ，若该枚弹药未失效，则具有截尾寿命时间 L =5 ），Y 和L 是独立的连续随机变量。设试验中 能且仅能观测到弹药样本的寿命I i = mi n｛Yi ，Li ｝， （i =1 ，…，7 ）. 若Yi ＞Li ，则第i 枚弹药在Li 后失 去观察（被截尾）。定义指示变量"：当Yi ＜Li ，"i = 1 ，当Yi ＞Li ，"i =0 . 于是得到弹药贮存寿命在随机 截尾试验中的似然函数为： L（"）= 7 （ ！ 7 - 】 7 i =1 "i ）！ H 7 i =1 （"）"［i 1 -F（"）］1 -"i . （4 ） 对于文献［1 ］中所指出的、贮存寿命服从指数分 布的弹药，设其贮存寿命"=（I 1 ，…，I 7 ）T 服从失效 率为#的指数分布，记为$（#），则其失效时间密度 （I i】#）= #exp（-I i #），可靠度及其失效时间分布 函数为：F（I i】#）= 1 -exp（-I i #）= 1 - R（I i】#）. 综上可得随机截尾条件下，关于$（#）的似然函数如 下： L（#】7 ，"，!）。 H 7 i =1 （I i】#）"i R（I i】#） （1 -"i ） 。 H 7 i =1 ［#iexp （-I i #i ）］" ［i exp （-I i #i ）］ （1 -"i ） . （5 ） 文献［1 ］中将贮存可靠度表示为故有可靠度的 条件贮存可靠度的乘积，认为固有可靠度应由弹药 出厂验收的试验数据得到。实际中，由于贮存环境 的复杂性使得该方法的局限性极强。由此引入失效 率#的伴随变量，建立指数回归模型，以确定贮存 寿命与某些环境变量（如：湿度、温度、辐射等）之间 的关系：令#i = exp（o T i"），式中：#i 为伴随变量，代 316 兵 工 学 报 第28 卷
第3期 基于MCMC稳态模拟方法的弹药贮存可靠性评估模型 317 表了影响弹药贮存寿命的主要因素：B为回归系数 向量，刻画了影响程度的大小。则指数回归模型的 (v(m-Dlog-exp(. (8) 似然函数如下： 同上，引入参数入的伴随变量建立Weibull回 L(BIn,t,x,)o I f(t;la)"R(o 归模型：令入；=exp(xB),则Weibull回归模型的似 然函数为 Ⅱ[exp(xB)exp(-t,exp(xB)]阝× L(m7lt,X,vcm防ep{[rB+ [exp(-t:exp(xB))]1-c v:(m -1)log(t;)-exp(x B)t] (9) ep(空xiBexp{-ew(B)}. (6) 3.2参数的贝叶斯估计 2.2参数的贝叶斯估计 由贝叶斯理论知，当m己知时，exp(入)的共轭 按照通常的处理方式，设B的先验分布为p维 先验密度服从Gamma分布：当m与入均未知时，尚 多元正态分布π（引40，S0),其中均值为40，协方差 没有合适的联合先验密度。这里令m服从Gamma 矩阵为S0,记为N,0,S0),则B后验分布为： 先验分布S(0,K0)、入服从正态先验分布V(u0, π（Bln,t,X,v)cL(Bln,t,X,v)π（Blo,So)c 品)，则m与入的联合后验分布如下： ep(空p)ep-会4eoxp小x π(m,入|n,t,X,D)cc L(m,aln,t,X,v)x(mlao,Ko)(loB)= cxp{←-B-nrs(B-小ac L(,λ|，t,X,v) e即（空B-会axB)- 1 aexp-(wo) √J2元 266 2(B-Ho)Sa(B-Ho). (7) 不难看出，(7)式的形式复杂，一般情况下，难以 i=l 给出其后验估计的精确形式，因而，对阝后验分布 的计算可借助前文中所介绍的MCMC模拟方法，利 e）-0m-2aA-a2. (10) 用Gibbs抽样方法从B的后验分布种抽样得出参数 引入参数入的伴随变量，得到Weibull回归模型m 的全条件分布密度： 与？的联合后验分布为 π（月|n,t,X,UB-)oc L(月，B-1n,t,X,v)×π，B-10S0). m,Bl,X,Em*em(空(gp+ 3 Veibull回归模型的贝叶斯分析 u,(m-1)log(t;）-exp(x:B)t）-kom- 3.1随机截尾的指数回归模型 2B-rs1(B-}. (11) 对于文献[1]中所指出的另一类贮存寿命服从 由上式也可以看出，π(m,入|n,t,X,v)的计算 形状参数为m、尺度参数为？的Weibull分布的弹 复杂，一般情况下难以给出后验结果精确形式，同上 药（记为W(m,7),记其失效时间密度为f(t,, 对后验分布的计算也可借助MCMC模拟方法。 7)=m-lexp(-t),可靠度R(tlm,7)及其 失效时间分布函数为： 4数值仿真分析 F(tilm,n)=1-exp(n )=1-R(tilm,n). 4.1仿真数据 令1=log(n),则上式化简为 假定某种弹药的贮存寿命服从指数分布e(入)， f(t;|m,a)=mtW-lexp(λ-exp(λ)t) 为了考察该种弹药在环境1与环境2中贮存寿命的 及：F(t;|m,A)=1-exp(-exp(入)t).则在随机 差异（即伴随变量一环境对其贮存寿命的影响，且设 截尾条件下，关于W(m,刀)的似然函数如下： 仅有两种环境可能)，分别取贮存于该两种环境的弹 L(m,aln,t,v)cc 药进行随机截尾试验（如2.1小节所述），并记录试 Ⅱflm,a)R(tlm,a)1-c 验中的射击数、失效数，弹药贮存后射击试验的仿真 数据如表1所示。设B=（B,B1)T,其中为截 m×n(以+ 距，31为环境变量x:的系数：设B~N2(0,10,I), 由此即得出如(6)式所示的一元指数回归模型
表了影响弹药贮存寿命的主要因素；!为回归系数 向量，刻画了影响程度的大小。则指数回归模型的 似然函数如下： L（! 7 ，!，"，"） 7 i =1 f（t i ） i R（t i ）（1 - i ） 7 i =1 ［exp （# T i!）exp （-t iexp （# T i!））］ i > ［exp （-t iexp （# T i!））］（1 - i ） exp｛ 7 i =1 i# T i!｝exp｛- 7 i =1 t iexp （# T i!）｝. （6 ） 2. 2 参数的贝叶斯估计 按照通常的处理方式，设!的先验分布为I 维 多元正态分布 （! #0 ，$0 ），其中均值为#0 ，协方差 矩阵为$0 ，记为 NI （#0 ，$0 ），则!后验分布为： （! 7 ，!，"，"） L（! 7 ，!，"，"） （! #0 ，$0 ） exp｛ 7 i =1 i# T i!｝exp｛- 7 i =1 t iexp （# T i!）｝> exp｛- 1 2 （!-#0 ）T $-1 0 （!-#0 ）｝ exp｛ 7 i =1 i# T i!- 7 i =1 t iexp （# T i!）- 1 2 （!-#0 ）T $-1 0 （!-#0 ）｝. （7 ） 不难看出，（7 ）式的形式复杂，一般情况下，难以 给出其后验估计的精确形式，因而，对!后验分布 的计算可借助前文中所介绍的 MCMC 模拟方法，利 用Gi bbs 抽样方法从!的后验分布种抽样得出参数 的全条件分布密度： （ j 7 ，!，"，"，! （-j ）） L（ j ，! （-j ） 7 ，!，"，"）> （ i ，! （-j ） #0 ，$0 ）. 3 Weibull 回归模型的贝叶斯分析 3. 1 随机截尾的指数回归模型 对于文献［1 ］中所指出的另一类贮存寿命服从 形状参数为 m 、尺度参数为 的 Wei bull 分布的弹 药（记为 W（m ， ）），记其失效时间密度为f（t i m ， ）= m t m -1 i exp（- t m i ），可靠度 R（t i m ， ）及其 失效时间分布函数为： F（t i m ， ）=1 -exp （ t m i ）=1 - R（t i m ， ）. 令 =log （ ），则上式化简为 f（t i m ， ）= mt m -1 i exp （ -exp （ ）t m i ） 及：F（t i m ， ）= 1 -exp （-exp （ ）t m i ）. 则在随机 截尾条件下，关于 W（m ， ）的似然函数如下： L（m ， 7 ，!，"） 7 i =1 f（t i m ， ） i R（t i m ， ） （1 - i ） m 7 i =1 i > exp｛ 7 i =1 i + 7 i =1 （ i （m -1 ）log （t i ）-exp （ ）t i m ｝. （8 ） 同上，引入参数 的伴随变量建立 Wei bull 回 归模型：令 i =exp （# T i!），则 Wei bull 回归模型的似 然函数为 L（m， 7 ，!，"，"） m 7 i =1 i >exp｛ 7 i =1 ［ i # T i!+ i （m -1 ）log （t i ）-exp （# T i!）t m i ］｝. （9 ） 3. 2 参数的贝叶斯估计 由贝叶斯理论知，当 m 已知时，exp（ ）的共轭 先验密度服从Ga mma 分布；当 m 与 均未知时，尚 没有合适的联合先验密度。这里令 m 服从 Ga mma 先验分布 （0 ， 0 ）、 服从正态先验分布N（ 0 ， 2 0 ），则 m 与 的联合后验分布如下： （m ， 7 ，!，"，"） L（m ， 7 ，!，"，"） （m 0 ， 0 ） （ 0 ， 2 0 ）= L（m ， 7 ，!，"，"）（ 0 0 （ 0 ）m 0 - exp （- 0 m ））> 1 2! 0 exp（- 1 2 2 0 （ - 0 ）2） m 0+ 7 i =1 i -1 exp｛ 7 i =1 i + 7 i =1 （ i （m -1 ）log （t i ）- exp （ ）t m i ）- 0 m - 1 2 2 0 （ - 0 ）2｝. （10 ） 引入参数 的伴随变量，得到 Wei bull 回归模型 m 与 的联合后验分布为 !（m ，! 7 ，!，"，"） m 0+ 7 i =1 i -1 exp｛ 7 i =1 （ i# T i!+ i （m -1 ）log （t i ）-exp （# T i!）t m i ）- 0 m - 1 2 （!-#0 ）T $-1 0 （!-#0 ）｝. （11 ） 由上式也可以看出， （m ， 7 ，t ，X ， ）的计算 复杂，一般情况下难以给出后验结果精确形式，同上 对后验分布的计算也可借助 MCMC 模拟方法。 4 数值仿真分析 4. 1 仿真数据 假定某种弹药的贮存寿命服从指数分布 （ ）， 为了考察该种弹药在环境1 与环境2 中贮存寿命的 差异（即伴随变量—环境对其贮存寿命的影响，且设 仅有两种环境可能），分别取贮存于该两种环境的弹 药进行随机截尾试验（如2. 1 小节所述），并记录试 验中的射击数、失效数，弹药贮存后射击试验的仿真 数据如表1 所示。设!=（ 0 ， 1 ）T，其中 0 为截 距， 1 为环境变量Ii 的系数；设!"N 2 （0 ，104 ，I ）， 由此即得出如（6 ）式所示的一元指数回归模型。 第3 期 基于 MCMC 稳态模拟方法的弹药贮存可靠性评估模型 317
318 兵工学报 第28卷 表1某种弹药贮存后射击试验的仿真数据 寿命差别不大：将环境变量x;及所求得的、31分 Tab.1 Simulated data in firing-test of certain 别代入模型，得到不同环境条件下弹药在各年份的 stored ammunition 贮存可靠度及可靠贮存寿命tR可表示为 贮存年 环境1 环境2 R(t:)=exp[-exp(-3.16+0.6125x;)t:], 限/a 射击数/次失效数/次 射击数/次失效数/次 InR(t;) 3 9 9 2 tR=-exp(-3.16+0.6125x) 5 5 2 5结束语 3 0 3 6 文献[1]将有效处理删失数据的生存分析方法 > 5 引入弹药贮存寿命的可靠性试验，使评估结果更接 近实际：而本文利用贝叶斯方法及回归的相关理论 4.2基于BUGS的仿真分析 大大缩小了试验样本的容量，通过对参数建立环境 BUGS是英国剑桥公共卫生研究所推出的利用 变量的回归模型，提高了模型的有效性，利用基于 MCMC方法进行贝叶斯推断的专用软件包，使用 Gibbs抽样的MCMC模拟方法与BUGS软件解决 BUGS可以很方便地对许多常用模型和分布进行 了该模型中高维数值计算的不便，提高了计算的精 Gibbs抽样。Gibbs抽样收敛后，可以得到参数的后 度，有利于该模型在可靠性分析理论中的推广。本 验分布的均值、标准差、95%置信区间等信息（见表 文仿真分析中的一元指数回归模型可以推广到多元 2),并给出后验分布的核密度估计图、参数的Gibbs 的情况中去，进而处理更为复杂的情况。需要指出 抽样动态图（见图1）等，使抽样结果更直观、可靠。 的是，先验分布的确定应视实际情况而定。不难看 这里基于表1所示的仿真数据，构建如上一元指数 出，该模型在弹药贮存寿命的可靠性分析中具有较 回归模型，截取前1000次迭代结果，从第1001次 单变量指数模型更为广泛的应用前景。 开始进行4000次迭代，部分显示结果如下。 参考文献(References) 表24000次抽样迭代的参数后验估计统计量 [1]李长福，夏建中，黄跃德，等.考虑删失数据时弹药贮存可靠性 Tab.2 Posterior estimate statistics of monitored 评估方法研究[].兵工学报，1996,17(4)：303-307. parameters with 4000 iterations LI Chang-fu,XIA Jian-zhong,HUANG Yue-de,et al.Method 参数 均值 标准差 2.5% 97.5% of assessment for the reliability of ammunition storage considering 号 -3.16 0.4337 -4.046 -2.403 censoring data[J].Acta Armamentarii,199,17(4):303-307. in Chinese) A 0.6125 0.5582 -0.5156 1.697 [2]刘金梅，王建萍，张力.生存分析法在弹药贮存可靠性评估中的 应用[」].弹箭与制导学报，2004,24(4)：332-334. LIU Jin-mei.WANG Jian-ping,ZHANG Li.The application of 件 survival analysis to evaluation of ammunition storage reliability [J].Joumal of Projectiles,Rockets,Missiles and Guidance, 1000 2000 3000 4000 5000 2004,24(4):332-334.(in Chinese). 迭代次数 [3] 郑波，张国安，李明.自然贮存环境下弹药系统贮存可靠性评估 [J门.质量与可靠性，2003，(3)：22-25 ZHENG Bo,ZHANG Guo-an,LI Ming.The research of ammu- nition storage reliability in natural environment.[.Quality and 广000 2000 3000 4000 5000 Reliability,2003.(3):22-25.(in Chinese) 达代次数 [4]郑波，李明，杨宝强，等.阵地环境下引信贮存可靠性研究[J门 探测与控制学报，2003,25(3)：51-53. 图14000次抽样选代后参数迭代的马尔可夫链轨迹 ZHENG Bo,LI Ming,YANG Bao-qiang,et al.The research of Fig.1 Trace plots of Markov chain for monitored fuze storage reliability in position environment.[J].Joumal of parameters after 4000 iterations Detection Control,2003,25(3):51-53.(in Chinese) [5]LawlessJ F,寿命数据中的统计模型与方法[M].茆诗松，译 从表2中可以看出：3的均值为-3.16,95% 北京：中国统计出版社，1998：37-38. 置信区间为(-4.046，-2.403)：B1的均值为 Lawless J F.Statistical models and methods for lifetime data[M]. 0.6125,95%置信区间为(-0.5156,1.697)，此区 MAO Shi-song,translated.Beijing:China Statistics Press, 间包含0点，可以认为在该两种环境下贮存的弹药 1998:37-38.(in Chinese)
表1 某种弹药贮存后射击试验的仿真数据 tab .1 Si mulated data i n firi ng-test of certai n stored a mmunition 贮存年 限／a 环境1 环境2 射击数／次 失效数／次 射击数／次 失效数／次 3 9 1 9 2 4 5 1 5 2 5 3 0 3 1 6 3 1 3 1 7 5 2 5 3 4. 2 基于BUGS 的仿真分析 BUGS 是英国剑桥公共卫生研究所推出的利用 MCMC 方法进行贝叶斯推断的专用软件包，使用 BUGS 可以很方便地对许多常用模型和分布进行 Gi bbs 抽样。Gi bbs 抽样收敛后，可以得到参数的后 验分布的均值、标准差、95% 置信区间等信息（见表 2 ），并给出后验分布的核密度估计图、参数的 Gi bbs 抽样动态图（见图1 ）等，使抽样结果更直观、可靠。 这里基于表1 所示的仿真数据，构建如上一元指数 回归模型，截取前1 000 次迭代结果，从第1 001 次 开始进行4 000 次迭代，部分显示结果如下。 表2 4 000 次抽样迭代的参数后验估计统计量 tab .2 Posterior esti mate statistics of monitored para meters wit h 4 000 iterations 参数 均值 标准差 2. 5% 97. 5% !0 -3. 16 0. 433 7 -4. 046 -2. 403 !1 0. 612 5 0. 558 2 -0. 515 6 1. 697 图1 4 000 次抽样迭代后参数迭代的马尔可夫链轨迹 Fig .1 trace plots of Markov chai n f or monitored para meters after 4 000 iterations 从表2 中可以看出：!0 的均值为-3. 16 ，95% 置 信 区 间 为（- 4. 046 ，- 2. 403 ）；!1 的 均 值 为 0. 612 5 ，95%置信区间为（-0. 515 6 ，1. 697 ），此区 间包含0 点，可以认为在该两种环境下贮存的弹药 寿命差别不大；将环境变量xi及所求得的!0 、!1 分 别代入模型，得到不同环境条件下弹药在各年份的 贮存可靠度及可靠贮存寿命rR 可表示为 r（ri）=exp ［-exp （-3. 16 + 0. 612 5xi）ri］， rR = - l nr（ri） exp （-3. 16 + 0. 612 5xi）. 5 结束语 文献［1 ］将有效处理删失数据的生存分析方法 引入弹药贮存寿命的可靠性试验，使评估结果更接 近实际；而本文利用贝叶斯方法及回归的相关理论 大大缩小了试验样本的容量，通过对参数建立环境 变量的回归模型，提高了模型的有效性，利用基于 Gi bbs 抽样的 MCMC 模拟方法与 BUGS 软件解决 了该模型中高维数值计算的不便，提高了计算的精 度，有利于该模型在可靠性分析理论中的推广。本 文仿真分析中的一元指数回归模型可以推广到多元 的情况中去，进而处理更为复杂的情况。需要指出 的是，先验分布的确定应视实际情况而定。不难看 出，该模型在弹药贮存寿命的可靠性分析中具有较 单变量指数模型更为广泛的应用前景。 参考文献（References ） ［1 ］ 李长福，夏建中，黄跃德，等. 考虑删失数据时弹药贮存可靠性 评估方法研究［J ］. 兵工学报，1996 ，17（4 ）：303 -307 . LI Chang-f u ，XI A Jian-zhong ，HUANG Yue- de ，et al . Met hod of assessment f or t he reliability of ammunition storage consi deri ng censori ng data ［J ］. Acta Ar mamentarii ，1996 ，17（4 ）：303 -307 . （i n Chi nese ） ［2 ］ 刘金梅，王建萍，张力. 生存分析法在弹药贮存可靠性评估中的 应用［J ］. 弹箭与制导学报，2004 ，24（4 ）：332 -334 . LI UJi n- mei ，WANG Jian- pi ng ，ZHANG Li . the application of survival analysis to evaluation of ammunition storage reliability ［J ］. Journal of Proj ectiles ，Rockets ，Missiles and Gui dance ， 2004 ，24（4 ）：332 -334 .（i n Chi nese ） ［3 ］ 郑波，张国安，李明. 自然贮存环境下弹药系统贮存可靠性评估 ［J ］. 质量与可靠性，2003 ，（3 ）：22 -25 . ZHENG Bo ，ZHANG Guo-an ，LI Mi ng . the research of ammunition storage reliability i n natural environment . ［J ］. @uality and Reliability ，2003 ，（3 ）：22 -25 .（i n Chi nese ） ［4 ］ 郑波，李明，杨宝强，等. 阵地环境下引信贮存可靠性研究［J ］. 探测与控制学报，2003 ，25（3 ）：51 -53 . ZHENG Bo ，LI Mi ng ，YANG Bao-Ciang ，et al . the research of f uze storage reliability i n position environment . ［J ］. Journal of Detection & Control ，2003 ，25（3 ）：51 -53 .（i n Chi nese ） ［5 ］ Lawless J F. 寿命数据中的统计模型与方法［M］. 茆诗松，译. 北京：中国统计出版社，1998 ：37 -38 . Lawless J F. Statistical models and met hods f or lif eti me data ［M］. MAO Shi-song ，translated . Beiji ng ：Chi na Statistics Press ， 1998 ：37 -38 .（i n Chi nese ） 318 兵 工 学 报 第28 卷
