孟慶東，楊 珍

(佛山市南海區(qū)技師學(xué)院，廣東 528237)

0 引言

對機(jī)電產(chǎn)品進(jìn)行可靠性分析最基本的工作是要確定其壽命分布參數(shù)。機(jī)電產(chǎn)品壽命服從兩參數(shù)Weibull 分布，其分布參數(shù)的貝葉斯估計(jì)需要對后驗(yàn)分布進(jìn)行二重積分。根據(jù)貝葉斯公式，在求解可靠性特征量MTBF 值區(qū)間估計(jì)時(shí)，需要進(jìn)行內(nèi)積分限為函數(shù)的二重積分。在進(jìn)行積分時(shí)，由于Bayes 后驗(yàn)分布復(fù)雜且內(nèi)積分限為函數(shù)，數(shù)值積分方法失效［1-3］。

利用Bayes 理論進(jìn)行統(tǒng)計(jì)推斷，需要積分運(yùn)算。對于兩參數(shù)Weibull 分布，需要對Bayes 后驗(yàn)分布進(jìn)行二重積分。一般有兩種方法進(jìn)行積分運(yùn)算，即數(shù)值積分與Monte Carlo 積分。

運(yùn)用數(shù)值積分計(jì)算時(shí)會(huì)遇到一些問題。當(dāng)被積函數(shù)較復(fù)雜時(shí)，將導(dǎo)致計(jì)算精度降低。同時(shí)，積分區(qū)間的大小也會(huì)直接影響計(jì)算結(jié)果，使得Bayes 后驗(yàn)計(jì)算出現(xiàn)不穩(wěn)定的趨勢。對于多重積分，因?yàn)槊總€(gè)外層積分都取決于內(nèi)層積分在一組點(diǎn)上的取值，所以會(huì)出現(xiàn)誤差的累積。

Monte Carlo 積分從目標(biāo)分布中抽樣模擬，利用隨機(jī)抽樣所得點(diǎn)進(jìn)行統(tǒng)計(jì)推斷，避免了數(shù)值積分中的求積節(jié)點(diǎn)計(jì)算與系統(tǒng)網(wǎng)格的劃分，使得積分計(jì)算的效率與可操作性顯著提高。但是，在Bayes 推斷中，運(yùn)用Monte Carlo 方法是有困難的，具體表現(xiàn)在Bayes 后驗(yàn)分布多為復(fù)雜、高維非標(biāo)準(zhǔn)的分布。當(dāng)目標(biāo)密度函數(shù)易于抽樣時(shí)，可以利用Monte Carlo 方法產(chǎn)生近似樣本，并利用樣本的函數(shù)值計(jì)算Bayes 后驗(yàn)期望估計(jì)。當(dāng)目標(biāo)密度函數(shù)難以抽樣時(shí)，可以利用MCMC 方法思想直接構(gòu)造出Bayes 后驗(yàn)分布π(m，η|T)的馬爾科夫鏈，基于所得樣本可以進(jìn)行各種統(tǒng)計(jì)推斷［4-6］。

1 MCMC 方法

MCMC 方法通過建立一個(gè)以后驗(yàn)分布π(x)為平穩(wěn)分布的馬爾科夫鏈，對π(x)進(jìn)行抽樣，然后根據(jù)所得樣本進(jìn)行統(tǒng)計(jì)推斷［7］。如通過抽樣得到π(x)的樣本X(1)，…，X(n)，則:

其中D 為狀態(tài)空間，于是可得估計(jì):

如果X(1)，…，X(n)是平穩(wěn)分布為π(x)的馬爾科夫過程的樣本時(shí)，上式也成立。

為了構(gòu)造馬爾科夫鏈，設(shè)給定狀態(tài)Xt條件下，下一步狀態(tài)為Xt+1。則任意選擇一個(gè)不可約轉(zhuǎn)移概率q(* ，* )以及一個(gè)函數(shù)α(* ，* )，定義p(x，y)=q(x，y)α(x，y)，則p(x，y)為馬氏鏈一個(gè)轉(zhuǎn)移核。如果鏈在時(shí)刻t 處于狀態(tài)x，即X(t)= x，首先由q(* | x)產(chǎn)生一個(gè)潛在的轉(zhuǎn)移x →y，然后根據(jù)概率α(x，y)決定是否轉(zhuǎn)移?？梢詮腢nif(0，1)中抽一個(gè)隨機(jī)數(shù)u，則:

通常，稱q(* ，* )為建議分布。

為使后驗(yàn)分布π(x)成為平穩(wěn)分布，在有q(* ，* )后，應(yīng)選擇一個(gè)α(* ，* )使相應(yīng)的p(x，y)以π(x)為其馬氏鏈平穩(wěn)分布，可以選擇:

此時(shí)，可以證明產(chǎn)生的馬爾科夫鏈?zhǔn)强赡娴模瑵M足平穩(wěn)方程，即: π(x)q(x，y)= π(y)q(y，x)

2 Weibull 分布下MTBF 值區(qū)間估計(jì)

2.1 先驗(yàn)分布為伽馬－ 逆伽馬分布的MTBF 值區(qū)間

選取均勻分布分布為建議分布，則:

其中，α(x，y)為Metropolis-Hastings 比率，my和ηy是在一步迭代中，從建議分布中抽取的建議值，mx和ηx是一步迭代的初始值。迭代過程如下:

①θ = (m，η)= (θ(t－1)0，θ(t－1)1);

②從建議分布Unif(－0.5，0.5)中產(chǎn)生候選點(diǎn);θ'

③令θ″ = θ + θ'，即讓建議值隨機(jī)游動(dòng)一個(gè)距離;

④根據(jù)上式計(jì)算接受概率α(θ，θ″);

⑤以概率α(θ，θ″)接受θ″ = θ(t)，否則令θ =θ(t)。

以某機(jī)電新產(chǎn)品運(yùn)行時(shí)間為現(xiàn)場樣本，用上述方法進(jìn)行計(jì)算，得到結(jié)果如表1 所示。

表1 MCMC 方法計(jì)算結(jié)果

圖1 和圖2 分別是參數(shù)m、η 以及MTBF 值的全部迭代軌跡和后驗(yàn)分布的密度直方圖。

圖1 參數(shù)m 和η 及MTBF 值迭代軌跡圖

圖2 參數(shù)m 和η 及MTBF 值的后驗(yàn)密度直方圖

由于計(jì)算可靠性特征量MTBF 值區(qū)間估計(jì)時(shí)，數(shù)值積分方法失效，依據(jù)工程經(jīng)驗(yàn)用參數(shù)點(diǎn)估計(jì)的區(qū)間估計(jì)端點(diǎn)來近似MTBF 值的置信區(qū)間。但是這種方法存在一定問題，一是求得的置信區(qū)間精度不高，有擴(kuò)大的趨勢;二是近似方法并沒有可靠的理論支撐。通過表1 的計(jì)算結(jié)果，可以看出MCMC 方法解決了上述問題，使得可靠性評(píng)估在質(zhì)量與效率上都得到提高。

2.2 先驗(yàn)分布為二元正態(tài)分布的MTBF 值區(qū)間估計(jì)

先驗(yàn)分布為二元正態(tài)分布時(shí)，計(jì)算結(jié)果如表2 所示。

表2 MCMC 方法計(jì)算結(jié)果

圖3 和圖4 分別是參數(shù)m、η 以及MTBF 值的全部迭代軌跡和后驗(yàn)分布的密度直方圖。

圖3 參數(shù)m 和η 及MTBF 值的迭代軌跡圖

圖4 參數(shù)m 和η 及MTBF 值的后驗(yàn)密度直方圖

3 結(jié)論

通過MCMC 構(gòu)建馬爾科夫鏈，使Bayes 后驗(yàn)計(jì)算運(yùn)用的數(shù)值積分方法轉(zhuǎn)化成從簡單的分布中抽樣并推斷?；诔闃铀梅植紖?shù)樣本，MTBF 值、可靠度、失效率等各種可靠性特征量的求解都可以有效實(shí)施，不再受數(shù)值積分的限制，提高了模型計(jì)算的穩(wěn)定性、可操作性與適應(yīng)性。MCMC 方法解決了數(shù)值積分問題帶來的不便，有利于貝葉斯理論在可靠性評(píng)估中的推廣。

［1］林靜. 基于MCMC 的貝葉斯生存分析理論及其在可靠性評(píng)估中的應(yīng)用［D］. 南京:南京理工大學(xué)，2008.

［2］李斌全，戴怡. 基于馬氏鏈蒙特卡洛方法的數(shù)控系統(tǒng)可靠性評(píng)估［J］. 組合機(jī)床與自動(dòng)化加工技術(shù)，2011(10):69－71.

［3］潘海濤，溫小霓. 基于MCMC 方法的GARCH 模型參數(shù)估計(jì)［J］. 統(tǒng)計(jì)與信息論壇，2009，24(4):12－16.

［4］曹小群，宋君強(qiáng)，張衛(wèi)民，等. 基于MCMC 方法的Lorenz混沌系統(tǒng)的參數(shù)估計(jì)［J］. 國防科技大學(xué)學(xué)報(bào)，2010，32(2):68－72，145.

［5］李凌，徐偉. 基于MCMC 方法的繼電器加速壽命試驗(yàn)分析［J］. 低壓電器，2010(2):13－16，35.

［6］王進(jìn)玲，曾聲奎，馬紀(jì)明，等. 混合變量系統(tǒng)基于MCMC的自適應(yīng)重要抽樣法［J］. 宇航學(xué)報(bào)，2012，33(1):94－101.

［7］黎光明，張敏強(qiáng). 先驗(yàn)信息對MCMC 方法估計(jì)概化理論方差分量變異量的影響［J］. 統(tǒng)計(jì)與決策，2012(7):27－30.