凌 光 戴 怡 李 曦

1.天津工程師范學(xué)院,天津,300222 2.華中科技大學(xué),武漢,430074

0 引言

Bayes方法是解決小子樣問題的一個重要方法,其核心和關(guān)鍵是如何利用先驗信息確定先驗分布。根據(jù)早期研究,數(shù)控系統(tǒng)失效時間服從雙參數(shù)Weibull分布。目前,雙參數(shù)先驗分布的確定未能得到很好解決：基于1/(θβ)形式的假設(shè)的方法[1]雖然簡化了雙參數(shù)相互間的關(guān)系,但在一定程度上影響了該方法的精度;獨立參數(shù)分布的截尾Γ分布和截尾Γ-1分布方法[2]計算復(fù)雜、不易操作,將對原分布的統(tǒng)計推斷的難度轉(zhuǎn)移到了參數(shù)分布計算上。

最大熵方法因充分利用了各種先驗信息并盡量避免了引入其他不確定因素而受到重視,許多學(xué)者為之作出不懈努力。Mazzuchi等[3]將非參數(shù)的 Dirichlet函數(shù)與最大熵相結(jié)合,提出了MED先驗分布;Kazama等[4]將最大熵方法應(yīng)用到不等式約束規(guī)劃問題中。Agrawal等[5]提出一種基于矩約束的最大熵雙參數(shù)聯(lián)合分布的確定方法,但該研究局限于一階矩約束。目前,關(guān)于一般形式下雙參數(shù)的最大熵先驗信息解的研究還鮮見報道。本文基于雙參數(shù)二階矩約束,求取了最大熵先驗信息解析解,并將之應(yīng)用于數(shù)控系統(tǒng)可靠性評估中。

1 利用最大熵方法確定聯(lián)合先驗分布

如上所述,數(shù)控系統(tǒng)失效時間服從雙參數(shù)Weibull分布 ,即

f(t;η,m)=m exp(-(t/η)m)/[η(t/η)m-1] t ＞0

它有η、m兩個參數(shù),其中,η為特征壽命(分布的0.632分位數(shù));m為形狀參數(shù),對密度函數(shù)形狀有很大影響。

根據(jù)二元聯(lián)合熵的定義和最大熵原理,確定雙參數(shù)聯(lián)合先驗分布 π(η,m),滿足

參數(shù)滿足的約束條件如下：

式中,ηi、mj分別為兩參數(shù)η、m的 i階和 j階原點矩;k、n分別為兩參數(shù)原點矩的最高階數(shù),可以相同,也可以不同。

對聯(lián)合密度函數(shù)π(η,m)的求解是個泛函意義下的條件極值問題。構(gòu)造如下輔助泛函：

式中,λ0、λ1、…、λk、c1、…、cn為拉格朗日乘子。

本文主要討論二階矩約束下雙參數(shù)聯(lián)合先驗分布 ,對于一般的k 、n,系數(shù)λ0、λ1、…、λk、c1、…、cn可通過數(shù)值計算[6]求得。故在二階矩約束下,即當(dāng)k=n=2時,先驗函數(shù)形式為

2 雙參數(shù)二階矩先驗分布的確定和穩(wěn)健性分析

利用最大熵方法求解雙參數(shù)聯(lián)合先驗分布過程中,需要參數(shù)樣本的二階矩信息。由于試驗所得數(shù)據(jù)為數(shù)控系統(tǒng)失效時間,故本文利用自助法,通過構(gòu)造再生子樣來獲取兩參數(shù)η、m的二階矩估計。詳細過程可參見文獻[7]。

在應(yīng)用先驗分布前,必須分析其對后驗統(tǒng)計推斷結(jié)果的影響,即分析先驗分布的穩(wěn)健性。一種簡單常用的先驗分布穩(wěn)健性判斷方法[8-9]就是利用相對似然邊際分布函數(shù)m*(t|π)。該相對似然函數(shù)邊際分布值代表以π(η,m)為先驗分布時現(xiàn)場樣本出現(xiàn)的相對概率。當(dāng)獲得現(xiàn)場樣本(t1,t2,…,tN)后,可計算m*(ti|π)(i=1,…,N)。如果所得值都不是非常小(一般不小于10-3),則說明先驗分布符合現(xiàn)場數(shù)據(jù)條件,即可以認為先驗分布π(η,m)是穩(wěn)健的。

對雙參數(shù)Weibull分布,聯(lián)合先驗分布 π(η,m)的相對似然邊際函數(shù)為

3 最大熵先驗信息解在數(shù)控系統(tǒng)可靠性評估中的應(yīng)用

據(jù)Bayes原理,設(shè)有數(shù)控系統(tǒng)壽命數(shù)據(jù)T=(t1,t2,…,tN),則其似然函數(shù)可表示為

利用已知先驗形式,則參數(shù)η、m聯(lián)合后驗分布可表示為

借助于數(shù)值計算,可以得到后驗分布的點估計：

最后根據(jù)Weibull分布性質(zhì)可知,數(shù)控系統(tǒng)平均后驗無故障工作時間的估計為

設(shè)有數(shù)控系統(tǒng)失效數(shù)據(jù)30個,如表1所示,利用這30個歷史數(shù)據(jù),應(yīng)用自助法生成200個容量為1000的自助樣本組,即求得η和m的參數(shù)序列,進一步可以計算兩參數(shù)的數(shù)學(xué)期望和方差：

表1 數(shù)控系統(tǒng)壽命數(shù)據(jù) h

η=3811.6,m=1.4552,σ2η=7920.0,σ2m=0.0011

從而該Weibull分布尺度參數(shù)與位置參數(shù)的雙參數(shù)聯(lián)合先驗密度函數(shù)為

設(shè)有現(xiàn)場數(shù)控系統(tǒng)的壽命數(shù)據(jù)20個,如表2所示。利用該數(shù)據(jù)對先驗分布的穩(wěn)健性進行驗證,計算出其相對似然邊際函數(shù)值亦如表2所示。從表2中數(shù)值可以看出,該先驗分布對于所得的現(xiàn)場數(shù)據(jù)而言,穩(wěn)健性良好,即利用該分布作為先驗分布是可取的。將現(xiàn)場數(shù)據(jù)形成的似然函數(shù)乘以由歷史數(shù)據(jù)求解得到的先驗分布,得到由式(8)計算的后驗分布,并可求出其點估計：

故而數(shù)控系統(tǒng)平均后驗無故障工作時間的估計為

表2 現(xiàn)場數(shù)控系統(tǒng)壽命數(shù)據(jù)及其相對似然函數(shù)值

4 結(jié)語

最大熵方法能夠最大限度地利用歷史樣本信息,適用于產(chǎn)品壽命長、實驗數(shù)據(jù)難以獲得的高可靠性產(chǎn)品的評估。本文提出的基于二階矩約束、通過最大熵原理確定的先驗分布,經(jīng)檢驗具有良好的穩(wěn)健性。算例表明,該方法適用于數(shù)控系統(tǒng)可靠性評估。

[1] Ferdous J,Uddin M U,Pandey M.Reliability Estimation with Weibull Inter Failure Times[J].Reliability Engineering and System Safety,1995,50：285-296.

[2] 張湘平.小子樣統(tǒng)計推斷與融合理論在武器系統(tǒng)評估中的應(yīng)用研究[D].長沙：國防科學(xué)技術(shù)大學(xué),2003.

[3] Mazzuchi T A,Soofi E S,Soyer R.Computation of Maximum Entropy Dirichlet for Modeling Lifetime Data[J].Computational Statistics&Data Analysis,2000,32：361-378.

[4] Kazama J,Tsujii J.Maximum Entropy Models with Inequality Constraints：a Case Study on Text Categorization[J].Machine Learning,2005,60：159-194.

[5] Agrawal D,Singh J K,Kumar A.Maximum Entropy-based Conditional Probability Distribution Runoff Model[J].Biosystems Engineering,2005,90(1)：103-113.

[6] 夏新濤,陳曉陽,張永振,等.航天軸承摩擦力矩的最大熵概率分布與bootstrap推斷[J].宇航學(xué)報,2007,28(5)：1395-1400.

[7] 張金槐,劉琦,馮靜.Bayes試驗分析方法[M].長沙：國防科技大學(xué)出版社,2007.

[8] Berger JO.統(tǒng)計決策論及貝葉斯分析[M].賈乃光,譯.北京：中國統(tǒng)計出版社,1998.

[9] 余禮軍,Ben Heydecker.基于概率加權(quán)矩與最大熵原理的交通流統(tǒng)計分布估計方法[J].公路交通科技,2008,25(2)：113-133.