基于改進兩階段法和自助法的壽命數(shù)據(jù)分析

王國東,牛占文,曲　亮,何　楨

(1.天津大學管理與經(jīng)濟學部,天津300072;2.航空經(jīng)濟發(fā)展河南省協(xié)同創(chuàng)新中心,河南鄭州450015)

基于改進兩階段法和自助法的壽命數(shù)據(jù)分析

王國東1,2,牛占文1,曲亮1,何楨1

(1.天津大學管理與經(jīng)濟學部,天津300072;2.航空經(jīng)濟發(fā)展河南省協(xié)同創(chuàng)新中心,河南鄭州450015)

加速壽命試驗的樣本量通常很小而且截尾會很嚴重,這樣會使得求得的極大似然統(tǒng)計量產(chǎn)生很大的偏差,進而影響到估計的準確度和精度.本文對兩階段法做了較大改進:證明了當各組數(shù)據(jù)服從Weibull分布且數(shù)據(jù)類型為無截尾或type II截尾時,存在兩個關(guān)于形狀參數(shù)和尺度參數(shù)的樞軸量,并采用無偏因子法修正極大似然統(tǒng)計量;考慮到異方差性,用加權(quán)最小二乘法替代最小二乘法估計加速模型的系數(shù).根據(jù)兩階段方法不能采用Fisher信息矩陣計算分位數(shù)置信區(qū)間的缺點,用自助法獲得置信區(qū)間.本文通過一個實例闡述改進的兩階段法的分析過程.另外,本文用相對偏差和均方誤差根作為評判分析方法的準則,將改進的兩階段法和極大似然法在不同的分位點處的壽命估計作了對比.仿真結(jié)果表明,改進的兩階段法在低分位點處的壽命估計更優(yōu).

加速壽命試驗;自助法;無偏因子;加權(quán)最小二乘法;Weibull分布

1　引　言

隨著產(chǎn)品性能的不斷提升,可靠性已經(jīng)成為質(zhì)量的一個重要維度.產(chǎn)品的可靠性水平直接決定了售后維修策略以及在保證期所花費的費用.生產(chǎn)商在確定產(chǎn)品價格前要定量評估產(chǎn)品的可靠性,然而在有限的時間內(nèi),采用一般的壽命試驗很難獲得充分的可靠性信息,一些產(chǎn)品在規(guī)定時間內(nèi)并沒有發(fā)生失效.因此,生產(chǎn)商一般會采用加速壽命試驗(accelerated life test,ALT)技術(shù)來獲得更多的可靠性信息.加速壽命試驗是在失效機理不變的基礎(chǔ)上,通過尋找產(chǎn)品壽命與應(yīng)力之間的數(shù)學關(guān)系-加速模型,利用加速應(yīng)力水平下的壽命特征去外推評估正常應(yīng)力下的壽命特征[1].

在有限的時間內(nèi),即使采用加速壽命試驗也不能獲得全部樣品的壽命信息.在高應(yīng)力水平下,安排的樣本量很小且不能保證全部失效;在一些較低的應(yīng)力水平下,樣品截尾可能會十分嚴重.截尾的樣品數(shù)據(jù)只是說明在這個時間點樣品沒有失效,并不知道確切的失效時間.另外,壽命數(shù)據(jù)一般不服從正態(tài)分布,而是服從一些偏態(tài)的分布,比如lognormal分布,gamma分布和Weibull分布.極大似然方法能夠很好的處理截尾和偏態(tài)數(shù)據(jù)[2-4],但是在小樣本或者樣品截尾很嚴重的情況下,極大似然估計的偏差會很大.Freeman[5]用仿真方法研究分組數(shù)據(jù)Weibull分布形狀參數(shù)極大似然估計的偏差.仿真結(jié)果表明,只有在每一組的樣本量達到20～40的時候形狀參數(shù)的極大似然估計才是近似無偏的.為了降低極大似然統(tǒng)計量的偏差,Yang等[6]用參數(shù)正交化的方法對分組數(shù)據(jù)情況下的極大似然統(tǒng)計量進行了修正.對于沒有截尾和type II截尾的數(shù)據(jù), Yang等[6]得到的統(tǒng)計量是樞軸量,這樣可以用無偏因子法進一步降低統(tǒng)計量的偏差.單組情況的無偏因子法已經(jīng)有很多學者進行了研究.Thoman等[7]證明在沒有截尾情況下Weibull分布形狀參數(shù)的極大似然統(tǒng)計量與形狀參數(shù)之比是樞軸量,并用Monte Carlo方法計算無偏因子值.Ross[8]給出了在沒有截尾數(shù)據(jù)情況下無偏因子的計算公式.Ross[9]給出了在type II截尾數(shù)據(jù)小于50%的情況下無偏因子的計算公式.然而,目前還沒有文獻對分組數(shù)據(jù)的無偏因子法進行研究.

錢萍等[10]提出用兩階段方法進行參數(shù)估計,即第一階段用極大似然法獲得分布參數(shù)的估計;第二階段用最小二乘法估計加速模型系數(shù).但是他們提出的兩階段法沒有考慮極大似然統(tǒng)計量的偏差,也沒有考慮到異方差性(heteroscedasticity).由于兩階段方法得到的參數(shù)估計并不是出自同一個似然函數(shù),故不能用Fisher信息矩陣方法來推斷置信區(qū)間.錢萍等[10]采用自助法(bootstrap method)獲得置信區(qū)間.

本文對兩階段法進行了改進,通過一個實例闡述如何用改進的兩階段法進行數(shù)據(jù)分析,并通過仿真的方法比較兩階段法和極大似然法在不同分位點處的壽命估計結(jié)果.

2　假　設(shè)

本文考慮包括m個應(yīng)力水平的加速壽命試驗X=(X1,X2,...,Xm)T,其中第i組的應(yīng)力水平是Xi.從一批產(chǎn)品中隨機抽取N個樣品,并分為m組樣本,各組樣本的樣本量分別為n1,n2,...,nm,且n1+ n2+···+nm=N.將第i組樣本放在應(yīng)力水平Xi下進行試驗.本文用到的假設(shè)有

1)每個應(yīng)力水平下樣品壽命都服從Weibull分布,其概率密度函數(shù)和分布函數(shù)分別為

其中ηi＞0,βi＞0分別是第i組的尺度參數(shù)和形狀參數(shù).tij＞0是第i組第j個樣品的失效時間.

2)尺度參數(shù)ηi的大小依賴于應(yīng)力水平Xi,加速模型為ln(ηi)=Xiγ,其中γ為加速模型的系數(shù)向量.

3)不同應(yīng)力水平下Weibull分布的形狀參數(shù)值相同,即β1=β2=···=βm=β.

3　改進的兩階段法

3.1改進兩階段法的步驟

階段1用極大似然法獲得參數(shù)估計

令tij表示第i組第j個樣品的壽命,則m組的聯(lián)合似然函數(shù)可以表示為[2]

如果樣品失效則δij=1;如果樣品沒有失效則δij=0.C是與數(shù)據(jù)截尾類型有關(guān)的常數(shù).

對數(shù)似然函數(shù)為

似然方程組可以表示為

參數(shù)β和ηi的估計值可以通過最大化式(2)得到,也可以通過求解似然方程組(3)得到.

階段2用加權(quán)最小二乘法獲得加速模型系數(shù)的估計

由式(2)可以得出關(guān)于參數(shù)β和ηi的二階偏導(dǎo)的負數(shù)為

另外,關(guān)于ηi和ηj所有的二階偏導(dǎo)均為0.

因此,觀測Fisher信息矩陣可以寫成

通過對觀測Fisher信息矩陣取逆可以求得極大似然統(tǒng)計量的漸進協(xié)方差矩陣,從而可以對應(yīng)得出i方差的估計值.這里用作為ln(i)方差的估計.

加速模型的參數(shù)估計公式為

其中

3.2用無偏因子法降低統(tǒng)計量偏差

極大似然統(tǒng)計量是漸進無偏的.當樣本量很大時,模型參數(shù)統(tǒng)計量的偏差是可以忽略的,但是如果樣本量很小或者截尾很嚴重,則必須考慮修正統(tǒng)計量的偏差.本小節(jié)介紹用無偏因子法降低統(tǒng)計量的偏差.

無偏因子法是依靠無偏因子來降低統(tǒng)計量的偏差.無偏因子是通過樞軸量獲得的.Thoman等[7]證明對于未分組情況下,當數(shù)據(jù)類型為沒有截尾或type II截尾時/β和β(log()-log(η))是樞軸量,即它的分布與參數(shù)β和η的值無關(guān),只是依賴于樣本量n和樣品失效的個數(shù)r.本文將證明這個結(jié)論對于多應(yīng)力水平的分組數(shù)據(jù)同樣成立.

3.3分位數(shù)估計和置信區(qū)間

生產(chǎn)商關(guān)心的是低分位點產(chǎn)品壽命的估計值和置信區(qū)間.由于產(chǎn)品的壽命很長,通常不會在使用條件下進行試驗,而是采用高于使用條件下的多個應(yīng)力水平進行試驗,并外推得到使用條件下產(chǎn)品的壽命特征信息.產(chǎn)品在使用條件下的壽命分布參數(shù)是通過加速模型估計得到的,而兩階段法并不能用統(tǒng)一的似然函數(shù)來獲得這些參數(shù)估計.因此采用Fisher信息矩陣的方法獲得分位數(shù)的置信區(qū)間并不可行.

Efron等[11]指出自助法(bootstrap method)是一種重抽樣方法,用來估計很難用解析方法計算得到的統(tǒng)計性質(zhì).下面介紹用自助法獲得置信區(qū)間的步驟并考慮到對統(tǒng)計量進行修偏.

4　實例

Zelen[12]設(shè)計了一個兩因子的部分壽命試驗來研究溫度應(yīng)力和電應(yīng)力對電容器壽命的影響,其中溫度應(yīng)力有2個水平,電應(yīng)力有4個水平,試驗分為8組,每組有8個樣品,試驗終止條件采用type II截尾,當某組有4個樣品失效時終止這組試驗.試驗的設(shè)計和壽命數(shù)據(jù)見表1.

表1　試驗設(shè)計及壽命數(shù)據(jù)Table 1 Test design and lifetime data

4.1分析

現(xiàn)在許多成熟的商業(yè)軟件可以對加速壽命試驗數(shù)據(jù)進行分析,常用的軟件有Minitab,R和SAS[13].這些軟件都是用極大似然法進行參數(shù)估計.本小節(jié)將分析比較極大似然法和兩階段法的結(jié)果.

采用改進兩階段法對電容器壽命數(shù)據(jù)進行分析,仿真次數(shù)為10 000次.兩階段法的第一階段分析結(jié)果見表2.

表2　階段一分析結(jié)果Table 2 Stage one analysis results

兩階段法的第二階段分析結(jié)果見表3;采用極大似然方法的分析結(jié)果見表4.

表3　階段二分析結(jié)果Table 3 Stage two analysis results

表4　極大似然法分析結(jié)果Table 4 Maximum likelihood method analysis results

對比表3和表4,不難看出對于加速模型的系數(shù)估計差別并不大,但極大似然估計的標準誤差更小.由于標準誤與β成反比[2],而極大似然法過大的估計了β,這使得表4中標準誤差的值過小.另外,表4中各系數(shù)估計的置信區(qū)間長度也較短.這一方面是由于標準誤差值較小;另一方面是由于得到的置信限是基于樣本量充分大的假設(shè),而這在本例中顯然并不可取.表3得到的置信限則不需要有關(guān)樣本量的假設(shè).

改進的兩階段法和極大似然法的分位點壽命估計和置信區(qū)間分別見表5和表6.本文選取的分位點p=0.01,0.05,0.10和0.50.從表5和表6看出,極大似然法在低分位點得到的估計值比兩階段法在低分位點的估計值要大.這主要是由于極大似然法沒有對統(tǒng)計量降偏造成的.極大似然法估計p分位點壽命的公式為,其中(p)表示標準最小極值分布分布函數(shù)的逆.當估計低分位點產(chǎn)品的壽命時,的值是負數(shù),而極大似然法過大的估計了β,這樣得到的p就會偏大.另外,極大似然法得到的置信區(qū)間也較短,這是因為極大似然法求解置信區(qū)間的一個假設(shè)是分位點處的壽命估計的標準差已知,顯然這在本試驗中并不滿足.如果樣本量很大,也可以認為標準差已知.但是在本試驗中每組只有8個樣品且有50%截尾.在這種情況下,如果使用極大似然法計算置信區(qū)間,會使得置信區(qū)間長度變短.表5中的置信區(qū)間結(jié)果是通過自助法獲得的.自助法不需要滿足上面提到的假設(shè),故結(jié)果更可信.

表5　改進的兩階段法的分位數(shù)參數(shù)估計和置信區(qū)間Table 5 Improved two-stage analysis percentile estimates and confidence intervals

表6　極大似然法的分位數(shù)參數(shù)估計和置信區(qū)間Table 6 Maximum likelihood analysis percentile estimates and confidence intervals

4.2仿真

仿真的結(jié)果見表7和表8,其中RB1和RB2分別表示極大似然方法和改進的兩階段方法的相對誤差(百分比);RMSE1和RMSE2分別表示極大似然方法和改進的兩階段方法的均方誤差根(百分比).從表7可以看出這兩種方法的相對誤差都隨著分位點的增大而減小,而改進的兩階段法的相對誤差在低分位點更小些;從表8可以看出均方誤差根隨著分位點的增大而增加,而兩階段法在低分位點的均方誤差根要更小些.從表7和表8可以看出,兩階段法在估計低分位點產(chǎn)品壽命時效果更好.

表7　分位數(shù)估計的相對偏差Table 7 Relative bias of the percentiles estimates

表8　分位數(shù)估計的均方誤差根Table 8 Root mean square error of the percentiles estimates

5　結(jié)束語

本文改進了加速壽命試驗可靠性數(shù)據(jù)分析的兩階段法:采用無偏因子法降低極大似然統(tǒng)計量的偏差,并用加權(quán)最小二乘法替代最小二乘法來估計加速模型的系數(shù).本文通過實例詳細闡述了統(tǒng)計量的偏差對產(chǎn)品分位點壽命的影響.文章將相對偏差和均方誤差根作為評判準則,并通過仿真方法來研究應(yīng)用兩階段法進行分位點壽命估計的性質(zhì),并與極大似然法進行比較.仿真結(jié)果表明,兩階段法在估計低分位點的壽命時要優(yōu)于極大似然法.而低分位點的壽命估計恰恰是工程技術(shù)人員最為關(guān)心的.

本文應(yīng)用仿真對不同分位點兩種方法壽命估計的結(jié)果進行了對比.但仿真的結(jié)果只是適用于特定的試驗條件(試驗包括8個應(yīng)力水平,每個應(yīng)力水平有8個樣品,采用type II截尾,當某組有4個樣品失效則終止這一組的試驗).對于樣品總量N以及每組樣本量ni和截尾比例對統(tǒng)計量偏差的影響本文沒有研究.這將是下一步要開展的工作.

企業(yè)的試驗設(shè)備都很有限,一般都會把同一個應(yīng)力水平的樣品安排在同一臺設(shè)備中進行試驗.這就會有子抽樣(subsampling)的問題[14].當存在子抽樣時,極大似然法將不再適用,但本文提出的兩階段法仍然適用.Freeman等[15]和Liu等[16]都提到兩階段法過大估計了形狀參數(shù).本文通過無偏因子法很好的解決了修偏問題.但是如果壽命試驗存在子抽樣問題,就不能用文中提到的自助法直接求得置信區(qū)間.如何改進自助法,使其能夠計算帶有子抽樣壽命數(shù)據(jù)分位數(shù)的置信區(qū)間,這也是下一步要開展的工作.

[1]姜同敏,王曉紅,袁宏杰,等.可靠性試驗技術(shù).北京:北京航空航天大學出版社,2012. Jiang T M,Wang X H,Yuan H J,et al.Reliability Test Techniques.Beijing:Beihang University Press,2012.(in Chinese)

[2]Meeker W Q,Escobar L A.Statistical Methods for Reliability Data.New York:John Wiley,1998.

[3]Lawless J F.Statistical Models and Methods for Lifetime Data.2nd Edition.New Jersey:John Wiley,2003.

[4]Nelson W.Accelerated Testing:Statistical Models,Test Plans and Data Analysis.New Jersey:John Wiley,2004.

[5]Freeman L J.A cautionary tale:Small sample size concerns for grouped lifetime data.Quality Engineering,2011,23(2):134–141.

[6]Yang Z,Lin D K J.Improved maximum-likelihood estimation for the common shape parameter of several Weibull populations. Applied Stochastic Models in Business and Industry,2007,23(5):373–453.

[7]Thoman D R,Bain L J,Antle C E.Inferences on the parameters of the Weibull distribution.Technometrics,1969,11(3):445–460.

[8]Ross R.Formulas to describe the bias and standard deviation of the ML-estimated Weibull shape parameter.IEEE Transactions on Dielectrics and Electrical Insulation,1994,1(2):247–253.

[9]Ross R.Bias and standard deviation due to Weibull parameter estimation for small data sets.IEEE Transactions on Dielectrics and Electrical Insulation,1996,3(1):28–42.

[10]錢萍,陳文華,李星軍,等.產(chǎn)品可靠性的Bootstrap回歸統(tǒng)計分析方法.儀器儀表學報,2010,31(11):2549–2554. Qian P,Chen W H,Li X J,et al.Bootstrap regression statistical analysis method for product reliability.Chinese Journal of Scientific Instrument,2010,31(11):2549–2554.(in Chinese)

[11]Efron B,Tibshirani R.An Introduction to the Bootstrap.New York:Chapman and Hall,1993.

[12]Zelen M.Factorial experiments in life testing.Technometrics,1959,1(3):269–288.

[13]Rigdon S E,Englert B R,Lawson I A,et al.Experiments for reliability achievement.Quality Engineering,2012,25(1):54–72.

[14]Freeman L J,Vining G G.Reliability data analysis for life test experiments with subsampling.Journal of Quality Technology,2010, 42(3):233–241.

[15]Freeman L J,Vining G G.Reliability data analysis for life test designed experiments with sub-sampling.Quality and Reliability Engineering International,2013,29(4):509–519.

[16]Liu X,Tang L C.Analysis of reliability experiments with blocking.Quality Technology and Quantitative Management,2013,10(2): 141–160.

Lifetime data analysis based on improved two-stage approach and bootstrap approach

Wang Guodong1,2,Niu Zhanwen1,Qu Liang1,He Zhen1
(1.College of Management and Economics,Tianjin University,Tianjin 300072,China 2.Cooperative Innovation Center for Avation Economy Development,Zhengzhou 450015,China)

When the sample sizes of accelerated life test are small and the data are heavily censored,the bias may be very serious and may affect the accuracy and precision of maximum likelihood estimates.In this paper,the two-stage approach is improved in two aspects:after proving that there are two pivotal quantities with shape parameter and scale parameters of Weibull distributions for complete data or type II censored data,Monte Carlo method is used to obtain the unbiased factors and reduce the bias of estimators;and least square method is replaced with the weighted least square method because of the heteroscedasticity.The two-stage approach cannot obtain confidence intervals of percentiles via Fisher information matrix.Therefore,the confidence intervals are predicted using bootstrap approach.After that,an example is provided to illustrate the procedures of proposed approach.Furthermore,the improved two-stage method is compared with maximum likelihood method based on relative bias and root mean square error criteria.The simulated results show that the improved two-stage method is better in the case of low percentiles.

accelerated life test;bootstrap method;unbiasing factor;weighted least square method;Weibull distribution

TB114.3

A

1000-5781(2016)05-0710-09

10.13383/j.cnki.jse.2016.05.015

2013-11-07;

2014-10-16.

國家自然科學基金資助項目(71402118;71071107);國家自然科學基金重點資助項目(70931004;71532008);國家自然科學基金杰出青年科學基金資助項目(71225006);航空科學基金資助項目(2013ZG55024;2014ZG55021);河南省高等學校重點科研資助項目(17A630072).

王國東(1981—),男,天津薊縣人,博士,講師,研究方向:試驗設(shè)計與可靠性改進,Email:gdwang@tju.edu.cn;

牛占文(1966—),男,內(nèi)蒙古赤峰人,博士,教授,研究方向:精益生產(chǎn)和質(zhì)量工程,Email:zw.niu@163.com;

曲亮(1982—),女,天津人,博士,講師,研究方向:質(zhì)量管理與可靠性分析,Email:qulucky@yahoo.com;

何楨(1967—),男,河南濮陽人,博士,教授,研究方向:質(zhì)量工程與六西格瑪管理,Email:zhhe0321@163.com.