劉平清

（六盤水師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，貴州 六盤水 553001）

壽命數(shù)據(jù)經(jīng)常出現(xiàn)在許多可靠性研究中，如加速壽命測(cè)試研究。可靠性是衡量產(chǎn)品壽命的一個(gè)重要指標(biāo)，指數(shù)分布族在可靠性分析中起著重要的作用，常用于分析壽命數(shù)據(jù)的分布有伯爾（Burr）-XII［1］分布、對(duì)數(shù)-浴缸（log-bathtub）分布［2］、韋布爾（Weibull）分布［3］和雙參數(shù)指數(shù)分布［4］，其中雙參數(shù)指數(shù)分布被廣泛應(yīng)用于壽命數(shù)據(jù)的建模，且不同的應(yīng)力模式影響試驗(yàn)的周期和精度。陳循等［5］對(duì)加速壽命試驗(yàn)做了較為系統(tǒng)的探討，指出步降應(yīng)力加速壽命試驗(yàn)方法應(yīng)用于與耗損型失效對(duì)應(yīng)的長(zhǎng)壽命評(píng)估可以大幅度提高試驗(yàn)效率，而不影響加速壽命試驗(yàn)的評(píng)估精度。英特哈布（Intekhab）和阿里夫（Arif）等［6］采用步降應(yīng)力加速壽命試驗(yàn)檢驗(yàn)機(jī)組壽命，利用漸進(jìn)式二類截尾的最大似然估計(jì)（Maximum Likelihood Estima-tion，MLEs），得到了漸近方差矩陣和協(xié)方差矩陣，構(gòu)造了估計(jì)量的置信區(qū)間，并進(jìn)行了仿真研究，驗(yàn)證了結(jié)果的準(zhǔn)確性。王錄和祁濤［7］基于步降應(yīng)力加速退化試驗(yàn)方法，對(duì)繞線聯(lián)軸器扭矩滯后的性能退化數(shù)據(jù)的收集及分析，采用Weibull壽命分布和逆冪率加速壽命模型，通過與恒應(yīng)力試驗(yàn)結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)兩種應(yīng)力模式下得到的壽命值基本相等，但前者大大縮短了試驗(yàn)周期。在加速壽命試驗(yàn)中，常用估計(jì)參數(shù)方法有極大似然估計(jì)（Maximum Likelihood Estimation，MLE）、貝葉斯（Bayes）方法、最優(yōu)線性無偏估計(jì)量（Best Linear Unbiased Evaluation，BLUE）和矩估計(jì)等，MLE在可靠性參數(shù)估計(jì)中有較廣的應(yīng)用，隨著計(jì)算機(jī)技術(shù)的增強(qiáng)，Bayes方法在可靠性試驗(yàn)中的應(yīng)用也得到了極大提升，特別是馬爾可夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）算法的引入，對(duì)Bayes方法在復(fù)雜問題中的應(yīng)用起到了極大的推動(dòng)作用。K.海霞（HAIXIA）和W.孔源（KONGYUAN）［8］利用逆高斯（Inverse Gaussian，IG）過程建立步降應(yīng)力加速退化試驗(yàn)（Accelerated Degradation Test，ADT）的評(píng)價(jià)方法及其可靠性函數(shù)，最大似然估計(jì)方法（MLE）和馬爾可夫鏈蒙特卡羅（MCMC）估計(jì)方法來估計(jì)所提出的退化模型中的未知參數(shù)。張鵬和劉平清［9］基于非參數(shù)統(tǒng)計(jì)符號(hào)檢驗(yàn)估計(jì)分位數(shù)區(qū)間值。榮格（Jung）和永（Yong）［10］為了克服漸進(jìn)式混合截尾方案的局限性，給出了一種魯棒Bayes方法來估計(jì)廣義Ⅰ型漸進(jìn)式混合截尾樣本的雙參數(shù)指數(shù)分布的未知參數(shù)。塔哈尼（Tahani）［11］在Ⅰ型混合截尾數(shù)據(jù)下威布爾指數(shù)分布（Weibull Exponential Distribution，WED）情況，利用MLE的漸近分布和M-H算法進(jìn)行參數(shù)的估計(jì)，推導(dǎo)出了基于平方誤差損失函數(shù)的貝葉斯估計(jì)。

綜上所述，現(xiàn)在研究中基于壽命數(shù)據(jù)的參數(shù)估計(jì)模型種類繁多，針對(duì)雙參數(shù)指數(shù)分布及步進(jìn)加速壽命試驗(yàn)和步降加速壽命試驗(yàn)的研究也不少，但是很少有學(xué)者在參數(shù)的估計(jì)中選擇MH算法對(duì)Bayes估計(jì)中的參數(shù)做估計(jì)。本文假定壽命分布服從雙參數(shù)指數(shù)分布，試驗(yàn)進(jìn)程采用連續(xù)觀測(cè)，而隨機(jī)移走發(fā)生在觀測(cè)到有產(chǎn)品失效時(shí)，且隨機(jī)移走數(shù)服從概率P一定的二項(xiàng)分布，結(jié)束機(jī)制為定數(shù)截尾，并利用Bayes估計(jì)方法對(duì)參數(shù)做估計(jì)。在估計(jì)過程中，采用MH算法，并以接受建議分布產(chǎn)生值的概率和失效數(shù)在不同應(yīng)力下的分配為優(yōu)化變量，以參數(shù)估計(jì)精度為優(yōu)化準(zhǔn)則，對(duì)試驗(yàn)進(jìn)行優(yōu)化設(shè)計(jì)。

1 模型建立

雙參數(shù)指數(shù)分布是壽命分布的常用分布之一，其分布函數(shù)和密度函數(shù)為：

其中，λ＞0為尺度參數(shù)，μ＞0為位置參數(shù)。

在壽命試驗(yàn)中，假設(shè)共有h+1個(gè)應(yīng)力水平，分別用Sh＞Sh-1＞...＞S0表示（其中S0為正常應(yīng)力水平），且有：Sh＞Sh-1＞...＞S0。并給定樣品總數(shù)為n，總失效產(chǎn)品數(shù)為m，在Si下的失效數(shù)為mi，即，在試驗(yàn)結(jié)束時(shí)，剩余未失效產(chǎn)品當(dāng)作移走。在Si下尺度參數(shù)和位置參數(shù)分別為λi和μi。用k i表示應(yīng)力Si下的加速因子，即，則應(yīng)力Si下的密度函數(shù)為：

將應(yīng)力Si+1轉(zhuǎn)換為Si的時(shí)刻為τh-i(i=0,1,2,...,h-1)，且定義τ0=0。設(shè)轉(zhuǎn)換為應(yīng)力S i時(shí)的試驗(yàn)產(chǎn)品數(shù)為ni，可繪制應(yīng)力Si下試驗(yàn)進(jìn)程，如表1所示：

2 步降加速壽命的統(tǒng)計(jì)分析

尼爾森（Nelson）在研究步進(jìn)壽命試驗(yàn)時(shí)提出了累積失效模型（Cumulative Exposure Model，CEM），通過模型的建立，可以得到該模型下的累積失效分布函數(shù)：

證明：用數(shù)學(xué)歸納法［12］證明如下：

證明1當(dāng)i=1時(shí)，已證推論成立；

證明2假設(shè)i=v成立，其中v∈[1,h)，即：

滿足歸納假設(shè)。

由證明1、證明2可知，對(duì)v∈[1,h)有推論成立。在Bayes估計(jì)中，先驗(yàn)分布的確定是一個(gè)非常重要的問題。常見的方法有：共軛先驗(yàn)分布、無信息先驗(yàn)分布、最大熵原則等［8］。假設(shè)已知位置參數(shù)μ0先驗(yàn)信息為E(μ0)=μ,var(μ0)=σ2，則可以考慮用最大熵原則確定μ0先驗(yàn)分布。而加速因子k i沒有任何先驗(yàn)信息，且服從貝葉斯假設(shè)。尺度參數(shù)λ0的共軛先驗(yàn)分布為伽瑪（Gamma）分布，且有，則各參數(shù)先驗(yàn)分布可確定為：

通過上文對(duì)模型的建立及說明，再根據(jù)貝葉斯理論可以得到模型的聯(lián)合先驗(yàn)密度函數(shù)：

并有后驗(yàn)密度函數(shù)：

由聯(lián)合先驗(yàn)密度函數(shù)形式的復(fù)雜決定了后驗(yàn)密度函數(shù)的復(fù)雜，通過經(jīng)驗(yàn)貝葉斯估計(jì)參數(shù)基本不可能實(shí)現(xiàn)，因此選用MH算法（MH算法是基于MCMC算法的優(yōu)化方法）對(duì)各參數(shù)進(jìn)行估計(jì)。則可以得到各個(gè)參數(shù)的滿條件分布的核分別為：

從各個(gè)參數(shù)的滿條件分布的核可以發(fā)現(xiàn)，λ0的滿條件分布為Gamma分布，μ0和加速因子ki的滿條件分布都不是常見的某個(gè)特定分布，并且相關(guān)函數(shù)形式的隨機(jī)數(shù)也不易于生成，針對(duì)本文模型，利用MH算法對(duì)參數(shù)進(jìn)行估計(jì)。關(guān)于MH算法的理論及設(shè)計(jì)過程［11］，MH算法是通過構(gòu)造馬爾可夫鏈生成貝葉斯后驗(yàn)分布樣本的典型算法，具體步驟如下：

步驟1初始化馬爾科夫鏈的狀態(tài)θi=θ1，令i=1,i=1,2,...,N,N為總的迭代次數(shù)；

步驟2從建議分布q中抽取一個(gè)候選樣本值θ*；

步驟3計(jì)算候選樣本值θ*的接受概率

步驟4從均勻分布U(0,1)中生成一個(gè)隨機(jī)數(shù)u，若u≤α(θi,θ*)，則接受候選樣本值θ*，令θi+1=θ*，否則令θi+1=θi；

步驟5令i=i+1，重復(fù)步驟2—步驟4，直到得到一個(gè)平穩(wěn)分布的馬爾科夫鏈。

在理論分析中，常需要對(duì)試驗(yàn)進(jìn)程進(jìn)行優(yōu)化，以期達(dá)到試驗(yàn)?zāi)茉谳^低成本下得到較為精確的結(jié)論的目的。在可靠性壽命試驗(yàn)中，常用的優(yōu)化變量有試驗(yàn)樣本的分配情況、不同應(yīng)力的時(shí)間轉(zhuǎn)換點(diǎn)等。對(duì)試驗(yàn)設(shè)計(jì)進(jìn)行優(yōu)化時(shí)，理論模型本身并沒有太大變化，但估計(jì)的精度和可信度卻能得到較大的提升。在優(yōu)化設(shè)計(jì)分析中，一般需要借助計(jì)算機(jī)軟件編程，給優(yōu)化變量在某個(gè)范圍內(nèi)一個(gè)較小的步長(zhǎng)，再分析不同取值下的估計(jì)情況，根據(jù)優(yōu)化準(zhǔn)則，找到優(yōu)化變量的最優(yōu)值，進(jìn)而實(shí)現(xiàn)對(duì)試驗(yàn)進(jìn)程的優(yōu)化。

3 案例分析

仿真目的是比較估計(jì)的效果。在試驗(yàn)中給定h=1，即只有兩個(gè)應(yīng)力水平，分別為S0,S1，將模型簡(jiǎn)化是為了對(duì)簡(jiǎn)單步降應(yīng)力模型的分析。而各參數(shù)的真實(shí)值為：μ0=10 000,λ0=0.001,k1=5。參加試驗(yàn)的產(chǎn)品總數(shù)N=60，截尾數(shù)n=20，并給定隨機(jī)移走概率為p=0.05，并且根據(jù)先驗(yàn)信息認(rèn)為μ=9 500，σ=500,α=10,β=10 000，用α1表示對(duì)建議分布生成值的接受率。由于本文考慮失效產(chǎn)品數(shù)在不同應(yīng)力水平下的分配問題，而不同分配下的壽命數(shù)據(jù)并不相同，鑒于數(shù)據(jù)量較大，在此只給出隨機(jī)移走數(shù)。隨機(jī)移走數(shù)據(jù)如表2所示：

表2 隨機(jī)移走數(shù)據(jù)

然后，利用R語言設(shè)計(jì)MH算法對(duì)模型進(jìn)行分析?？梢怨烙?jì)出參數(shù)在不同情況下的值，結(jié)果如表3所示。

數(shù)品產(chǎn)效失的到測(cè)觀下S1平水應(yīng)在18171615141312111098765432 10453.60 10467.67 10394.34 10471.15 10428.89 10476.71 10454.67 10426.25 10419.29 10491.19 10440.26 10482.18 10382.81 10431.46 10472.24 10450.84 10479.51 u^o 0.001544 0.001647 0.001700 0.001647 0.001690 0.001615 0.001722 0.001628 0.001701 0.001571 0.001706 0.001667 0.001736 0.001714 0.001746 0.001638 0.001667 λ^0 α 1=0.2 4.709477 4.70288 4.704793 4.701539 4.69946 4.713836 4.711166 4.699492 4.694968 4.695484 4.687406 4.706858 4.690946 4.692237 4.704553 4.701649 4.69659 k^1 10127.00 10132.88 10165.13 10109.01 10140.14 10163.60 10143.75 10167.97 10146.57 10116.37 10190.34 10105.20 10139.52 10158.52 10175.42 10107.59 10165.52 u^0 0.001405 0.001436 0.001418 0.001413 0.001425 0.001312 0.001335 0.001377 0.001359 0.001378 0.001429 0.001209 0.001328 0.001328 0.001331 0.001433 0.001440 λ^0 α 1=0.5 4.859740 4.862617 4.878955 4.882453 4.855749 4.845094 4.874400 4.879009 4.863204 4.873815 4.862143 4.857071 4.864703 4.869073 4.849664 4.853769 4.874169 k^1 9946.728 9960.941 9943.811 9945.332 9943.645 9955.955 9958.052 9959.054 9931.038 9936.703 9935.234 9935.043 9955.053 9941.955 9949.447 9943.806 9950.536 u^0 0.001112 0.001153 0.001126 0.001104 0.001121 0.001134 0.001107 0.001177 0.001142 0.001175 0.001148 0.001155 0.001157 0.001117 0.001120 0.001140 0.001150表3不同α1取值在S1下不同失效產(chǎn)品數(shù)的參數(shù)估計(jì)值α1取值待估參數(shù)λ^0 α 1=0.8 5.055736 5.063368 5.063282 5.076198 5.051946 5.056438 5.069496 5.041554 5.042923 5.032718 5.056483 5.085819 5.039721 5.056938 5.042168 5.054540 5.051472 k^1 9793.490 9777.956 9813.362 9799.009 9799.635 9786.812 9825.653 9797.054 9782.417 9818.236 9791.128 9803.985 9829.451 9797.716 9777.640 9793.925 9786.737 u^0 0.001002 0.000997 0.001004 0.001016 0.001015 0.001010 0.000998 0.000998 0.000993 0.000992 0.001001 0.001000 0.001009 0.000988 0.001013 0.000990 0.000993 λ^0 α 1=1 5.214507 5.199333 5.206135 5.183910 5.178126 5.202639 5.219798 5.215460 5.198949 5.201915 5.214155 5.203924 5.209938 5.207259 5.178622 5.221747 5.163864 k^1

可用圖形表示趨勢(shì)如圖1所示：

圖1 不同a1取值在S1下不同失效產(chǎn)品數(shù)的估計(jì)值趨勢(shì)

從表3可以看出，在對(duì)u0做估計(jì)時(shí)，當(dāng)a1一定時(shí)，隨著應(yīng)力水平S1下分配的失效樣品的增加，u0的估計(jì)值有一定的波動(dòng)，但波動(dòng)較小，而且沒有明顯的變化趨勢(shì)。然而a1變化時(shí)，u0的估計(jì)值的波動(dòng)較大，并且在a1=0.8時(shí)估計(jì)值達(dá)到最優(yōu)。對(duì)加速因子k i的估計(jì)值變化情況與u0有相同的結(jié)論。在對(duì)l0做估計(jì)時(shí)，當(dāng)a1一定時(shí)，隨著應(yīng)力水平S1下分配的失效樣品的增加，l0的估計(jì)值的波動(dòng)形式與u0一致，但是在a1=1時(shí)，l0的估計(jì)值最優(yōu)。通過分析可知，由于l0的滿條件分布是已知分布，因此a1越大估計(jì)結(jié)果越好，而對(duì)滿條件分布不是某特定分布的參數(shù)而言，在利用MH算法對(duì)參數(shù)做估計(jì)時(shí)，a1值并不是越大越好，也不是越小越好，具體分析時(shí)需要根據(jù)建議分布對(duì)滿條件分布的擬合程度進(jìn)行更深入的考慮。

4 結(jié)語

通過模型推導(dǎo)和模擬分析，在雙參數(shù)指數(shù)分布的步降壽命試驗(yàn)中，失效數(shù)在不同應(yīng)力水平下的分配對(duì)參數(shù)估計(jì)精度的影響較??；但是當(dāng)u0,λ0,k1,…,k h其中一個(gè)參數(shù)的滿條件分布不是特定分布，且該分布隨機(jī)數(shù)不易生成時(shí)，對(duì)建議分布生成值的接受率做合理選擇，能有效地提高參數(shù)估計(jì)的精度。通過對(duì)模型的分析，在Bayes估計(jì)中選擇MH算法對(duì)參數(shù)做估計(jì)，并利用案例數(shù)據(jù)分析簡(jiǎn)單步降應(yīng)力下不同a1（接受建議分布生成值的概率）和m1（應(yīng)力水平S1下失效樣本數(shù)）對(duì)各參數(shù)估計(jì)精度的影響，發(fā)現(xiàn)m1對(duì)參數(shù)的估計(jì)精度影響不大，但a1對(duì)參數(shù)的估計(jì)精度有一定的積極意義。提出的雙參數(shù)指數(shù)分布模型中的思想不僅適用步降應(yīng)力實(shí)驗(yàn)，同時(shí)也適用于壽命服從Weibull分布其他類型壽命實(shí)驗(yàn)的可靠性分析。