孟慶東,楊 珍
(佛山市南海區(qū)技師學院,廣東 528237)
對機電產(chǎn)品進行可靠性分析最基本的工作是要確定其壽命分布參數(shù)。機電產(chǎn)品壽命服從兩參數(shù)Weibull 分布,其分布參數(shù)的貝葉斯估計需要對后驗分布進行二重積分。根據(jù)貝葉斯公式,在求解可靠性特征量MTBF 值區(qū)間估計時,需要進行內(nèi)積分限為函數(shù)的二重積分。在進行積分時,由于Bayes 后驗分布復雜且內(nèi)積分限為函數(shù),數(shù)值積分方法失效[1-3]。
利用Bayes 理論進行統(tǒng)計推斷,需要積分運算。對于兩參數(shù)Weibull 分布,需要對Bayes 后驗分布進行二重積分。一般有兩種方法進行積分運算,即數(shù)值積分與Monte Carlo 積分。
運用數(shù)值積分計算時會遇到一些問題。當被積函數(shù)較復雜時,將導致計算精度降低。同時,積分區(qū)間的大小也會直接影響計算結果,使得Bayes 后驗計算出現(xiàn)不穩(wěn)定的趨勢。對于多重積分,因為每個外層積分都取決于內(nèi)層積分在一組點上的取值,所以會出現(xiàn)誤差的累積。
Monte Carlo 積分從目標分布中抽樣模擬,利用隨機抽樣所得點進行統(tǒng)計推斷,避免了數(shù)值積分中的求積節(jié)點計算與系統(tǒng)網(wǎng)格的劃分,使得積分計算的效率與可操作性顯著提高。但是,在Bayes 推斷中,運用Monte Carlo 方法是有困難的,具體表現(xiàn)在Bayes 后驗分布多為復雜、高維非標準的分布。當目標密度函數(shù)易于抽樣時,可以利用Monte Carlo 方法產(chǎn)生近似樣本,并利用樣本的函數(shù)值計算Bayes 后驗期望估計。當目標密度函數(shù)難以抽樣時,可以利用MCMC 方法思想直接構造出Bayes 后驗分布π(m,η|T)的馬爾科夫鏈,基于所得樣本可以進行各種統(tǒng)計推斷[4-6]。
MCMC 方法通過建立一個以后驗分布π(x)為平穩(wěn)分布的馬爾科夫鏈,對π(x)進行抽樣,然后根據(jù)所得樣本進行統(tǒng)計推斷[7]。如通過抽樣得到π(x)的樣本X(1),…,X(n),則:
其中D 為狀態(tài)空間,于是可得估計:
如果X(1),…,X(n)是平穩(wěn)分布為π(x)的馬爾科夫過程的樣本時,上式也成立。
為了構造馬爾科夫鏈,設給定狀態(tài)Xt條件下,下一步狀態(tài)為Xt+1。則任意選擇一個不可約轉移概率q(* ,* )以及一個函數(shù)α(* ,* ),定義p(x,y)=q(x,y)α(x,y),則p(x,y)為馬氏鏈一個轉移核。如果鏈在時刻t 處于狀態(tài)x,即X(t)= x,首先由q(* | x)產(chǎn)生一個潛在的轉移x →y,然后根據(jù)概率α(x,y)決定是否轉移。可以從Unif(0,1)中抽一個隨機數(shù)u,則:
通常,稱q(* ,* )為建議分布。
為使后驗分布π(x)成為平穩(wěn)分布,在有q(* ,* )后,應選擇一個α(* ,* )使相應的p(x,y)以π(x)為其馬氏鏈平穩(wěn)分布,可以選擇:
此時,可以證明產(chǎn)生的馬爾科夫鏈是可逆的,滿足平穩(wěn)方程,即: π(x)q(x,y)= π(y)q(y,x)
選取均勻分布分布為建議分布,則:
其中,α(x,y)為Metropolis-Hastings 比率,my和ηy是在一步迭代中,從建議分布中抽取的建議值,mx和ηx是一步迭代的初始值。迭代過程如下:
①θ = (m,η)= (θ(t-1)0,θ(t-1)1);
②從建議分布Unif(-0.5,0.5)中產(chǎn)生候選點;θ'
③令θ″ = θ + θ',即讓建議值隨機游動一個距離;
④根據(jù)上式計算接受概率α(θ,θ″);
⑤以概率α(θ,θ″)接受θ″ = θ(t),否則令θ =θ(t)。
以某機電新產(chǎn)品運行時間為現(xiàn)場樣本,用上述方法進行計算,得到結果如表1 所示。
表1 MCMC 方法計算結果
圖1 和圖2 分別是參數(shù)m、η 以及MTBF 值的全部迭代軌跡和后驗分布的密度直方圖。
圖1 參數(shù)m 和η 及MTBF 值迭代軌跡圖
圖2 參數(shù)m 和η 及MTBF 值的后驗密度直方圖
由于計算可靠性特征量MTBF 值區(qū)間估計時,數(shù)值積分方法失效,依據(jù)工程經(jīng)驗用參數(shù)點估計的區(qū)間估計端點來近似MTBF 值的置信區(qū)間。但是這種方法存在一定問題,一是求得的置信區(qū)間精度不高,有擴大的趨勢;二是近似方法并沒有可靠的理論支撐。通過表1 的計算結果,可以看出MCMC 方法解決了上述問題,使得可靠性評估在質(zhì)量與效率上都得到提高。
先驗分布為二元正態(tài)分布時,計算結果如表2 所示。
表2 MCMC 方法計算結果
圖3 和圖4 分別是參數(shù)m、η 以及MTBF 值的全部迭代軌跡和后驗分布的密度直方圖。
圖3 參數(shù)m 和η 及MTBF 值的迭代軌跡圖
圖4 參數(shù)m 和η 及MTBF 值的后驗密度直方圖
通過MCMC 構建馬爾科夫鏈,使Bayes 后驗計算運用的數(shù)值積分方法轉化成從簡單的分布中抽樣并推斷。基于抽樣所得分布參數(shù)樣本,MTBF 值、可靠度、失效率等各種可靠性特征量的求解都可以有效實施,不再受數(shù)值積分的限制,提高了模型計算的穩(wěn)定性、可操作性與適應性。MCMC 方法解決了數(shù)值積分問題帶來的不便,有利于貝葉斯理論在可靠性評估中的推廣。
[1]林靜. 基于MCMC 的貝葉斯生存分析理論及其在可靠性評估中的應用[D]. 南京:南京理工大學,2008.
[2]李斌全,戴怡. 基于馬氏鏈蒙特卡洛方法的數(shù)控系統(tǒng)可靠性評估[J]. 組合機床與自動化加工技術,2011(10):69-71.
[3]潘海濤,溫小霓. 基于MCMC 方法的GARCH 模型參數(shù)估計[J]. 統(tǒng)計與信息論壇,2009,24(4):12-16.
[4]曹小群,宋君強,張衛(wèi)民,等. 基于MCMC 方法的Lorenz混沌系統(tǒng)的參數(shù)估計[J]. 國防科技大學學報,2010,32(2):68-72,145.
[5]李凌,徐偉. 基于MCMC 方法的繼電器加速壽命試驗分析[J]. 低壓電器,2010(2):13-16,35.
[6]王進玲,曾聲奎,馬紀明,等. 混合變量系統(tǒng)基于MCMC的自適應重要抽樣法[J]. 宇航學報,2012,33(1):94-101.
[7]黎光明,張敏強. 先驗信息對MCMC 方法估計概化理論方差分量變異量的影響[J]. 統(tǒng)計與決策,2012(7):27-30.