基于改進兩階段法和自助法的壽命數據分析

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

(1.天津大學管理與經濟學部,天津300072;2.航空經濟發展河南省協同創新中心,河南鄭州450015)

基于改進兩階段法和自助法的壽命數據分析

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

(1.天津大學管理與經濟學部,天津300072;2.航空經濟發展河南省協同創新中心,河南鄭州450015)

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

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

1　引　言

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

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

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

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

2　假　設

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

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

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

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

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

3　改進的兩階段法

3.1改進兩階段法的步驟

階段1用極大似然法獲得參數估計

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

如果樣品失效則δij=1;如果樣品沒有失效則δij=0.C是與數據截尾類型有關的常數.

對數似然函數為

似然方程組可以表示為

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

階段2用加權最小二乘法獲得加速模型系數的估計

由式(2)可以得出關于參數β和ηi的二階偏導的負數為

另外,關于ηi和ηj所有的二階偏導均為0.

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

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

加速模型的參數估計公式為

其中

3.2用無偏因子法降低統計量偏差

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

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

3.3分位數估計和置信區間

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

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

4　實例

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

表1　試驗設計及壽命數據Table 1 Test design and lifetime data

4.1分析

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

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

表2　階段一分析結果Table 2 Stage one analysis results

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

表3　階段二分析結果Table 3 Stage two analysis results

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

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

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

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

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

4.2仿真

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

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

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

5　結束語

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

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

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

[1]姜同敏,王曉紅,袁宏杰,等.可靠性試驗技術.北京:北京航空航天大學出版社,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]錢萍,陳文華,李星軍,等.產品可靠性的Bootstrap回歸統計分析方法.儀器儀表學報,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—),男,天津薊縣人,博士,講師,研究方向:試驗設計與可靠性改進,Email:gdwang@tju.edu.cn;

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

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

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