基于馬氏鏈蒙特卡洛方法

1、基于馬氏鏈蒙特卡洛方法的數(shù)控系統(tǒng)可靠性評估李斌全，戴怡 （天津職業(yè)技術(shù)師范大學(xué)，機(jī)械工程學(xué)院，天津 300222）摘要：數(shù)控系統(tǒng)屬于高可靠性產(chǎn)品，應(yīng)用貝葉斯理論，已經(jīng)成為其可靠性評估的重要方法。由于數(shù)控系統(tǒng)壽命試驗(yàn)數(shù)據(jù)服從Weibull分布，導(dǎo)致后驗(yàn)分布出現(xiàn)數(shù)值大、高維復(fù)雜的情況，貝葉斯計(jì)算依靠的數(shù)值積分方法難以實(shí)施。根據(jù)馬氏鏈蒙特卡洛（MCMC）方法思想，建立后驗(yàn)平穩(wěn)分布的馬爾科夫鏈，對分布參數(shù)的貝葉斯估計(jì)進(jìn)行求解，解決了數(shù)值積分問題，保證可靠性評估的有效實(shí)施。通過與BUGS軟件結(jié)果進(jìn)行比較，表明提高了模型計(jì)算的穩(wěn)健性、有效性及精度。關(guān)鍵詞：數(shù)控系統(tǒng)；馬氏鏈蒙特卡洛（MCMC）方法；可靠性評

2、估；Weibull分布；貝葉斯理論中圖分類號：TH 17The reliability evaluation of NC system based on Markov chain Monte Carlo methodLi Bin-quan，Dai Yi(School of Mechanical Engineering, Tianjin University of Technology and Education, Tianjin 300222，China)Abstract: NC system was high-reliability product, so the application o

3、f Bayesian theory had become an important means of its reliability evaluation. The failure data of NC system obeyed Weibull distribution which had complex forms. So the posterior distribution became more complicated and the numerical integration which Bayesian computing depended on was not available

4、. The problem of numerical integration was solved by Markov chain Monte Carlo (MCMC) method. It ensures the effective implementation of reliability evaluation. Compared with BUGS package, the result of Bayesian estimation shows that it increases the robustness and effectiveness of the calculation. K

5、ey words: NC system；Markov chain Monte Carlo (MCMC) method；reliability evaluation；Weibull distribution；Bayesian theory0 引言可靠性工程中，對產(chǎn)品進(jìn)行可靠性研究或可靠性特征量預(yù)測，最基本的工作是要確定其壽命分布參數(shù)。據(jù)多位學(xué)者的工作，數(shù)控系統(tǒng)壽命服從兩參數(shù)Weibull分布1,2。其分布參數(shù)的貝葉斯估計(jì)需要對后驗(yàn)分布進(jìn)行二重積分。根據(jù)貝葉斯公式，后驗(yàn)分布為基金項(xiàng)目：國家自然科學(xué)基金（50875186）；國家數(shù)控機(jī)床科技重大專項(xiàng)（2009ZX04014-013）（1）其中，和為W

6、eibull分布兩參數(shù)，為樣本似然函數(shù)，為先驗(yàn)分布。則分布參數(shù)的貝葉斯估計(jì)為 （2） （3）由于數(shù)控系統(tǒng)壽命較長且Weibull分布形式復(fù)雜，使得被積函數(shù)數(shù)值巨大。在實(shí)際評估中，計(jì)算（2）、（3）式的內(nèi)層積分時(shí)就已經(jīng)出現(xiàn)浮點(diǎn)數(shù)溢出情況，導(dǎo)致依靠被積函數(shù)值計(jì)算的數(shù)值積分方法失效隨著馬氏鏈蒙特卡洛（MCMC）方法的不斷發(fā)展，貝葉斯后驗(yàn)分布計(jì)算有望得到徹底解決。通過MCMC方法，復(fù)雜的后驗(yàn)分布被直觀的模擬出來，這就使得參數(shù)求解運(yùn)用的數(shù)值積分方法轉(zhuǎn)化成從簡單的分布中抽樣并推斷?；诔闃铀梅植紖?shù)樣本，后驗(yàn)分布的統(tǒng)計(jì)計(jì)算簡單易行，各種可靠性特征量求解的可操作性顯著提高。目前貝葉斯計(jì)算多用國外的BUGS

7、軟件進(jìn)行3,4，但是該軟件未公布其核心算法，這給掌握算法實(shí)質(zhì)以及進(jìn)行算法改進(jìn)帶來困難。另一方面BUGS軟件采用的Gibbs抽樣是MCMC方法的特例，所以自主掌握MCMC方法，并應(yīng)用于可靠性評估有重要意義。本文利用MCMC方法對參數(shù)的貝葉斯估計(jì)進(jìn)行求解，并與BUGS軟件結(jié)果進(jìn)行比較，表明其適用于數(shù)控系統(tǒng)的可靠性評估。1 基于MCMC方法的貝葉斯可靠性1.1截尾樣本的Weibull似然函數(shù)構(gòu)建對數(shù)控系統(tǒng)進(jìn)行有替換定時(shí)截尾壽命試驗(yàn)，得到，其中，為完全失效數(shù)據(jù)樣本，為截尾數(shù)據(jù)樣本。兩參數(shù)Weibull分布概率密度函數(shù)為，生存函數(shù)為，其中為形狀參數(shù)，為尺度參數(shù)。則截尾樣本似然函數(shù)為 （4）1.2 后驗(yàn)分

8、布形式的確立取參數(shù)的先驗(yàn)分布為伽馬分布，其中為形狀參數(shù)，為尺度參數(shù)。參數(shù)的先驗(yàn)分布為逆伽馬分布，其中為形狀參數(shù)，為尺度參數(shù)。則兩參數(shù)的聯(lián)合先驗(yàn)分布為 （5）由貝葉斯公式得到后驗(yàn)分布： （6）由于后驗(yàn)分布數(shù)值巨大，高維復(fù)雜，貝葉斯計(jì)算依靠的數(shù)值積分方法難以實(shí)施。但是可以利用MCMC方法思想解決上述問題。2 后驗(yàn)平穩(wěn)分布的馬爾科夫鏈建立貝葉斯后驗(yàn)量的計(jì)算可歸結(jié)為關(guān)于后驗(yàn)分布的積分5。 MCMC方法思想是通過建立后驗(yàn)平穩(wěn)分布的馬爾科夫鏈6，并從中抽樣，來獲取后驗(yàn)分布樣本進(jìn)行統(tǒng)計(jì)推斷，這就不需要對被積函數(shù)進(jìn)行數(shù)值計(jì)算，從而保證求解過程的有效實(shí)施和穩(wěn)定性。 令，表示從易于抽樣的建議分布產(chǎn)生一個潛在的轉(zhuǎn)移

9、，然后根據(jù)概率決定是否轉(zhuǎn)移。如果鏈在時(shí)刻t處于狀態(tài)x，即，從中抽一個隨機(jī)數(shù)u， （7）則此時(shí)，為馬氏鏈一個轉(zhuǎn)移核。在選取后，根據(jù)（6）式，概率取為 （8）這時(shí)，轉(zhuǎn)移核以為馬氏鏈平穩(wěn)分布。選取均勻分布分布為建議分布，此時(shí)， （9）其中，為轉(zhuǎn)移概率，和是在一步迭代中，從建議分布中抽取的建議值，和是一步迭代的初始值。過程如下:（1）（2）從建議分布中產(chǎn)生候選點(diǎn)（3）令，即讓建議值隨機(jī)游動一個距離（4）根據(jù)式（9）計(jì)算接受概率（5）以概率接受，否則令通過上述過程，后驗(yàn)分布的馬爾科夫鏈被建立起來，并得到Weibull分布參數(shù)和的樣本，基于所得樣本就可以對分布參數(shù)進(jìn)行各種統(tǒng)計(jì)推斷。3 仿真分析通過有替換定

10、時(shí)截尾壽命試驗(yàn)，獲得數(shù)控系統(tǒng)壽命的現(xiàn)場樣本如表1所示，表1數(shù)控系統(tǒng)壽命現(xiàn)場樣本（單位：h）完全失效數(shù)據(jù)176256339455495117612601896236221342133322431853158151120134102438013604251截尾數(shù)據(jù)11369315232231245625883322102242131021250246614503210310250116011145141522332根據(jù)表1數(shù)據(jù)，用上述方法編程進(jìn)行計(jì)算，得到結(jié)果如表2表2MCMC方法計(jì)算結(jié)果建議分布接受率E(m)E()sd(m)sd()均勻分布0.511.3540990.171.8BUGS是進(jìn)行貝葉

11、斯推斷的軟件包，用該軟件處理表1數(shù)據(jù)，得到結(jié)果如表3表3應(yīng)用BUGS軟件計(jì)算結(jié)果nodemeansdstartsamplem1.320.1810015700041114.1100157000通過表2計(jì)算結(jié)果，表明MCMC方法的有效性與穩(wěn)定性，與表3進(jìn)行比較，可以看出其計(jì)算精度也有所提高。圖1、2分別是m和的全部迭代軌跡和后驗(yàn)分布的密度直方圖。（a）參數(shù)m（b）參數(shù)圖1 參數(shù)m和的迭代軌跡圖（a）參數(shù)m（b）參數(shù)圖2 后驗(yàn)分布的密度直方圖4 數(shù)控系統(tǒng)可靠性特征量的計(jì)算用建議分布為均勻分布時(shí)建立的馬氏鏈對數(shù)控系統(tǒng)可靠性進(jìn)行估計(jì)，Weibull分布概率密度函數(shù)采用傳統(tǒng)方法計(jì)算MTBF值的區(qū)間估計(jì)是復(fù)

12、雜而且難以實(shí)施的。利用MCMC方法可以有效地予以估計(jì)。圖3是MTBF值的迭代軌跡。圖3 MTBF值的迭代軌跡圖可以得到，MTBF值的90%置信區(qū)間為3666.76，3989.30 。5 結(jié)論（1）應(yīng)用貝葉斯理論進(jìn)行可靠性評估是一種有前途的方法。由于以往不能很好地解決高維數(shù)值積分問題，從而制約了該理論在實(shí)際中的應(yīng)用。MCMC方法解決了高維數(shù)值積分問題帶來的不便，有利于貝葉斯理論在可靠性評估中的推廣。（2）由于后驗(yàn)分布被模擬構(gòu)造出來，后驗(yàn)分布的各種統(tǒng)計(jì)推斷都可從簡單分布中抽樣進(jìn)行?；诔闃铀梅植紖?shù)樣本，可靠度、失效率等各種可靠性特征量的求解都可以有效實(shí)施，不再受數(shù)值積分的限制，提高了模型計(jì)算的

13、穩(wěn)定性、可操作性與適應(yīng)性參考文獻(xiàn)1 戴怡，周云飛，陳學(xué)東等. 加工中心故障分布規(guī)律及其研究方法J. 系統(tǒng)工程與電子技術(shù)，2004，26(3)：413-415.2 張海波，賈亞洲，周廣文. 數(shù)控系統(tǒng)故障間隔時(shí)間分布模型的研究J. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2005，37(2)：198-200.3 林靜. 基于MCMC的貝葉斯生存分析理論及其在可靠性評估中的應(yīng)用D. 南京：南京理工大學(xué)，2008.4 林靜，韓玉啟，朱慧明等. 基于MCMC穩(wěn)態(tài)模擬的Weibull回歸模型及其可靠性應(yīng)用J.系統(tǒng)仿真學(xué)報(bào)，2006，18(5) ：1161-1163.5 茆詩松，王靜龍，濮曉龍. 高等數(shù)理統(tǒng)計(jì)M.北京：高等教育出版社，1998.6 Geof