疲勞統(tǒng)計(jì)學(xué)智能化的進(jìn)一步研究

徐家進(jìn)，高鎮(zhèn)同

1.力中國際融資租賃有限公司,廣州 510620 2.北京航天航空大學(xué) 航空科學(xué)與工程學(xué)院,北京 100083

《疲勞統(tǒng)計(jì)學(xué)智能化中的高鎮(zhèn)同法》中指出，威布爾分布是一種全態(tài)分布，即不僅可描寫左偏、右偏的數(shù)據(jù)，在一定程度上也可描寫對稱的數(shù)據(jù)。在這個(gè)意義上比高斯分布用途更為廣泛，特別是在擬合結(jié)構(gòu)的疲勞壽命方面起到非常重要的作用。但因在確定威布爾分布3參數(shù)方面遇到了困難，因此人們往往將“最小壽命”歸零。但這種做法是不合適的，因?yàn)椤白钚勖钡奈锢硪饬x是可靠度為100%的壽命，應(yīng)該稱為“安全壽命”，這是非常重要的。為了確定威布爾分布的3個(gè)參數(shù)，以往采取圖解法和解析法，前者使用不便，且誤差比較大；后者則是要解3個(gè)聯(lián)立的超越方程組，盡管可用計(jì)算機(jī)來解，但有時(shí)候會出現(xiàn)解出的安全壽命這個(gè)參數(shù)會大于實(shí)際數(shù)據(jù)的最小值，即會存在這種不自冾的問題。文獻(xiàn)[1]中利用最小二乘法和Python特點(diǎn)提出了高鎮(zhèn)同法，比較好地解決了這個(gè)問題。這是疲勞統(tǒng)計(jì)學(xué)智能化的開始。而在本文則是更進(jìn)一步，只要給出(疲勞壽命)樣本數(shù)據(jù)通過計(jì)算機(jī)就可自動判斷該數(shù)據(jù)的總體服從的最佳分布(正態(tài)分布或威布爾分布)，并同時(shí)給出該分布的參數(shù)。同時(shí)智能化地給出每一個(gè)疲勞壽命的置信區(qū)間。恰恰是疲勞統(tǒng)計(jì)學(xué)智能化的一種成果。

盡管解析法有時(shí)候會出現(xiàn)不自冾的問題，但在實(shí)際應(yīng)用中并非一定會出現(xiàn)不自冾問題。此時(shí)用高鎮(zhèn)同法也會得到一組威布爾分布3個(gè)參數(shù)的估計(jì)值，而與解析法得到的一組估計(jì)的參數(shù)值是不相同的。于是產(chǎn)生了到底采取哪一組參數(shù)更加合適的問題。這就需要給出一個(gè)所謂“優(yōu)劣”的標(biāo)準(zhǔn)。注意到高鎮(zhèn)同法的本質(zhì)就是要選擇適當(dāng)?shù)陌踩珘勖鼌?shù)使得得到的威布爾分布與“理想”的可靠度的相關(guān)系數(shù)最大。故可用相關(guān)系數(shù)絕對值的大小作為衡量標(biāo)準(zhǔn)是恰當(dāng)?shù)摹５浑y看出相關(guān)系數(shù)可在線性擬合中起作用，遇到非線性擬合時(shí)就必須要用決定系數(shù)大小來作為標(biāo)準(zhǔn)。這個(gè)標(biāo)準(zhǔn)也可作為衡量各算法結(jié)果“優(yōu)劣”的標(biāo)準(zhǔn)。進(jìn)一步還需研究各參數(shù)的“置信區(qū)間”問題。但研究表明，同時(shí)給出這3個(gè)參數(shù)的置信區(qū)間是有困難的。因此，本文不僅給出了解決擬合“優(yōu)劣”的標(biāo)準(zhǔn)問題，同時(shí)也給出了威布爾分布在某種可靠度下每一點(diǎn)(即疲勞壽命)的置信區(qū)間的一個(gè)智能化解決方案，而這恰恰也是人們在實(shí)際工作中最需要的。

1 相關(guān)系數(shù)和決定系數(shù)

一般情況下此文說到相關(guān)系數(shù)都是指皮爾遜相關(guān)系數(shù)(Pearson Correlation Coefficient)，若、為2個(gè)隨機(jī)變量，則它們的相關(guān)系數(shù)為

=Cov(,)(Var()Var())12

(1)

式中：Var()、Var()分別為和的方差;Cov(,)是協(xié)方差[(-())(-())]，(·)為期望。

由柯西不等式：

(2)

式中：、分別表示和的第個(gè)分量。若將和看成為2個(gè)矢量，那么它們的所謂相關(guān)系數(shù)實(shí)際上就是這2個(gè)矢量的夾角余弦，這就是相關(guān)系數(shù)的幾何意義。為零表明2者成90°(正交)，因此完全不相關(guān)。反過來||越大則2者的夾角越小因此越相關(guān)，=+1稱之為完全正相關(guān),而=-1完全負(fù)相關(guān)。

相關(guān)系數(shù)在統(tǒng)計(jì)學(xué)中有著非常重要的地位，是衡量2個(gè)變量線性相關(guān)密切程度的指標(biāo)，與此有關(guān)的另外一個(gè)指標(biāo)就是決定系數(shù)。注意到數(shù)據(jù)由某曲線擬合的程度可定義為

1-SSE/SST=1-Var(-)Var()

(3)

=+

(4)

式中：、分別為斜率和截距，由此可得

=

(4a)

=-

(4b)

因

(5a)

但由式(4)和式(4b)可得

-=(-)

(5b)

故有

-=(-)-(-)

(-)

(6)

注意到式(4a)中的定義可知式(6)=0，

(7)

=SSR/SST

(8)

由式(5b),

=[()]=

注意，式(1)可改寫為=()12。

不難看出，的范圍是[0,1],而的范圍是[-1,1]。另外的統(tǒng)計(jì)意義是擬合程度，也可理解為變量對變量的解釋程度(回歸誤差和)的大小。由此可見，決定系數(shù)這個(gè)概念的意義更加廣泛，只不過在線性擬合的情況下

=

(9)

在非線性擬合情況下所謂相關(guān)系數(shù)就沒有意義，決定擬合程度只能用決定系數(shù)來判斷，且只能用式(3) 這個(gè)定義式。

事實(shí)上在比較各種分布擬合數(shù)據(jù)的過程中會發(fā)現(xiàn)不同的“參照系”會得到不同的決定系數(shù)。比如用高鎮(zhèn)同法為了求出威布爾分布的位置參數(shù)先要確定與數(shù)據(jù)“理想可靠度”最大的決定系數(shù)(因?yàn)榇藭r(shí)是線性擬合，因此仍然可用相關(guān)系數(shù))。然后再由最小二乘法定出威布爾分布的其他2個(gè)參數(shù)(形狀參數(shù)和尺度參數(shù))。但要與解析法及高斯分布比較究竟哪一個(gè)方法擬合得更好呢？就得將這3種方法得到的擬合可靠度曲線和理想可靠度來比較。此時(shí)就不再是線性擬合的情況下，就不能再使用相關(guān)系數(shù)而必須要用決定系數(shù)了。

2 以理想可靠度作為基準(zhǔn)的決定系數(shù)

所謂“理想可靠度”，實(shí)際上是指“可靠度估計(jì)量”：

=1-(+1)

(10)

式中：為第個(gè)疲勞壽命的可靠度;為疲勞壽命總個(gè)數(shù)。

為什么稱之為“理想”可靠度？這是因?yàn)槭?10)對于任何分布都是適用的，因此提供了一個(gè)參照標(biāo)準(zhǔn)。實(shí)際上高鎮(zhèn)同法的基礎(chǔ)也是建立在與理想可靠度的相關(guān)性之上的，即以此為標(biāo)準(zhǔn)來選擇位置參數(shù)。不過此時(shí)用的是對數(shù)空間中的相關(guān)系數(shù)，即ln(ln(-1/))與ln(-)的相關(guān)系數(shù)。而不是直接拿與威布爾分布的可靠度[exp(-(-)/)]的決定系數(shù)來比較(、、分別為威布爾分布的位置、尺度和形狀參數(shù))。

用文獻(xiàn)[2]第136頁的疲勞壽命數(shù)據(jù)來說明這個(gè)問題。利用Python代碼，運(yùn)行后可得

={350, 380, 400, 430, 450, 470, 480, 500, 520, 540, 550, 570, 600, 610, 630, 650, 670, 730, 770, 840}，表示疲勞壽命數(shù)組，該數(shù)組的長度LS= 20， 疲勞壽命均值= 557.0， 標(biāo)準(zhǔn)差= 132.15， 中值= 545.0, 偏態(tài)系數(shù)=0.408。

用高鎮(zhèn)同法結(jié)果：= 2.040，= 320.98，= 276.60,= 0.999 218，= 0.605，,、,分別為相應(yīng)最大相關(guān)系數(shù)的形狀參數(shù)、長度參數(shù)及安全壽命(位置參數(shù)) 和相關(guān)系數(shù)。與理想的可靠度的決定系數(shù)來比較：高鎮(zhèn)同法計(jì)算的威布爾分布的決定系數(shù)=0.998 24，高斯分布的決定系數(shù)=0.989 48，用解析法計(jì)算：=2.211,=312.32，=280.39，=0.503，用解析法計(jì)算的威布爾分布的決定系數(shù)=0.995 86。

由此可見，解析法和高鎮(zhèn)同法得到的結(jié)果與理想可靠度的決定系數(shù)結(jié)果相差無幾，這從它們得到的3個(gè)參數(shù)相差“不大”也是可預(yù)料得到的。這也說明數(shù)據(jù)擬合并沒有一個(gè)絕對正確的答案。高斯分布與理想可靠度相應(yīng)的決定系數(shù)也是非常接近1的，盡管比威布爾分布相應(yīng)的決定系數(shù)要小一些。同時(shí)可注意一下偏態(tài)系數(shù)，原始數(shù)據(jù)的偏態(tài)系數(shù)為0.408，而用高鎮(zhèn)同法和解析法得到的偏態(tài)系數(shù)分別為0.605和0.503。盡管解析法得到的偏態(tài)系數(shù)與原始的偏態(tài)系數(shù)更加接近一點(diǎn)，但決定系數(shù)卻小了一點(diǎn)，若以決定系數(shù)為標(biāo)準(zhǔn)則應(yīng)取高鎮(zhèn)同法的結(jié)果。

至于文獻(xiàn)[1]中例2已經(jīng)表明解析法無法計(jì)算出不自冾的結(jié)果，但高鎮(zhèn)同法是可解的，且高鎮(zhèn)同法得到的相關(guān)系數(shù)要比高斯分布要大，故可采用高鎮(zhèn)同法的結(jié)果,看起來高鎮(zhèn)同法的結(jié)果更佳。但并沒有給出嚴(yán)格的論證，故要對威布爾分布要做進(jìn)一步的研究。

3 利用數(shù)據(jù)分組對威布爾分布做進(jìn)一步的研究

如何才能對威布爾分布做進(jìn)一步研究呢？考慮到樣本容量比較大的時(shí)候往往會對其分組，其實(shí)通常做統(tǒng)計(jì)分析時(shí)也常常需要就數(shù)據(jù)分組。數(shù)據(jù)分組的好處是將“大數(shù)據(jù)”變成為好處理的“小數(shù)據(jù)”，又能保留原來“大數(shù)據(jù)”的“最主要”的特性。不過為此是要付出代價(jià)的，代價(jià)之一是分組是有講究的，是需要人為“干預(yù)”的。另一個(gè)代價(jià)則是無論怎樣分組多多少少都會“扭曲原意”， 盡管如此分組不失為一個(gè)統(tǒng)計(jì)分析的方法。比如在文獻(xiàn)[2]第253頁介紹卡方檢驗(yàn)時(shí)就將在同一應(yīng)力水平下測得的100個(gè)試件的疲勞壽命成為9組，再根據(jù)9組中不夠10次的“組”進(jìn)行必要的“合并”成為了6組，然后再進(jìn)行卡方檢驗(yàn)。結(jié)果表明這組數(shù)據(jù)是不符合高斯分布的。那么屬于什么分布呢？是不是符合威布爾分布呢？這就需要進(jìn)一步研究。其實(shí)這也是發(fā)現(xiàn)高鎮(zhèn)同法的一個(gè)契機(jī)，這一點(diǎn)在文獻(xiàn)[1]中已經(jīng)披露了。只不過當(dāng)時(shí)沒有將此數(shù)據(jù)進(jìn)行分組而是直接利用高鎮(zhèn)同法對其進(jìn)行計(jì)算。

為了更加直觀地感受不同方法不同分布對于數(shù)據(jù)的擬合程度，修改Python代碼后得如下結(jié)果：

={3.08, 3.26, 3.32, 3.48, 3.49, 3.56, 3.69, 3.7, 3.78, 3.79, 3.8, 3.87, 3.95, 4.07, 4.08, 4.1, 4.12, 4.2, 4.24, 4.25, 4.28, 4.31, 4.31, 4.36, 4.54, 4.58, 4.6, 4.62, 4.63, 4.65, 4.67, 4.67, 4.72, 4.73, 4.75, 4.77, 4.8, 4.82, 4.84, 4.9, 4.92, 4.93, 4.95, 4.96, 4.98, 4.99, 5.02, 5.03, 5.06, 5.08, 5.06, 5.1, 5.12, 5.15, 5.18, 5.2, 5.22, 5.38, 5.41, 5.46, 5.47, 5.53, 5.56, 5.6, 5.61, 5.63, 5.64, 5.65, 5.68, 5.69, 5.73, 5.82, 5.86, 5.91, 5.94, 5.95, 5.99, 6.04, 6.08, 6.13, 6.16, 6.19, 6.21, 6.26, 6.32, 6.33, 6.36, 6.41, 6.46, 6.81, 7.0, 7.35, 7.82, 7.88, 7.96, 8.31, 8.45, 8.47, 8.79, 9.87}，數(shù)據(jù)未分組結(jié)果：= 5.32，= 1.29，=5.07，= 1.021。數(shù)據(jù)分組結(jié)果：= 6.79，= 2.09，= 6.77，=0.043， 每組上限:3.84，4.68，5.52，6.35，7.19，8.03，8.87，9.87，每組含有的個(gè)數(shù):11，21，29，25，5，4，4，1，累積頻率:0.999，0.89，0.68，0.39，0.14，0.09，0.05，0.01。

對分組高鎮(zhèn)同法：=1.948，=2.73，=3.76,=0.999 247，=0.666，對分組用解析法：估計(jì)=3.253，估計(jì)=4.26，估計(jì)=1.50，=0.091。 與分組得到的實(shí)際積累頻率的決定系數(shù)：由高鎮(zhèn)同法得到的決定系數(shù)=0.994 82，由解析法得到的決定系數(shù)=0.955 84，由高斯分布得到的決定系數(shù)=0.957 46。不分組高鎮(zhèn)同法：=2.147，=2.87，=2.78,=0.993 763，=0.539。與不分組的可靠度的決定系數(shù)：由高鎮(zhèn)同法得到的決定系數(shù)=0.989 03，由高斯分布得到的決定系數(shù)=0.979 11。

由上述計(jì)算結(jié)果可得出如下幾個(gè)看法：

1) 原來用解析法得出最低壽命會大于實(shí)際數(shù)據(jù)中的開始3個(gè)數(shù)據(jù)，這顯然是不符合“安全壽命”這個(gè)概念的(即不自冾)。而現(xiàn)在將這100個(gè)數(shù)據(jù)分成為8組后，這個(gè)問題表面上是不存在了。但圖1表明用解析法得到的結(jié)果與高斯分布幾乎到了不可分辨的地步，而且明顯偏離了實(shí)際數(shù)據(jù)，故解析法和高斯分布這2個(gè)方法都不可取。

2) 從圖1及相關(guān)系數(shù)都可發(fā)現(xiàn)，對于分組數(shù)據(jù)高鎮(zhèn)同法得到的3參數(shù)威布爾分布的擬合曲線都比用解析法得到的3參數(shù)威布爾分布及高斯分布擬合得要好。盡管如此，由高鎮(zhèn)同法得出的安全壽命仍然明顯大于原始疲勞壽命數(shù)據(jù)中的最小值。因此在這個(gè)意義上分組并不是一個(gè)好主意。

圖1 數(shù)據(jù)分組的各種分布與實(shí)際可靠度擬合情況Fig.1 Fitting of various distributions and ideal reliability of data grouping

3) 再回到不分組的情況。從曲線擬合圖(圖2只給出了20個(gè)數(shù)據(jù)而不是100個(gè)數(shù)據(jù)點(diǎn)，主要原因是后者的圖太密，以致看不清)，用高鎮(zhèn)同法計(jì)算出來的3參數(shù)威布爾分布與高斯分布也是區(qū)別“不大”的，很難說哪一條更好。也只能從決定系數(shù)這個(gè)角度去看，此時(shí)可確定用高鎮(zhèn)同法計(jì)算出來的3參數(shù)威布爾分布比高斯分布是略勝一籌的。而且從偏態(tài)系數(shù)這個(gè)角度看，原數(shù)據(jù)偏態(tài)系數(shù)是相當(dāng)大的，故用高斯分布來擬合確實(shí)是不合適的。即偏態(tài)系數(shù)可作為判斷是否符合高斯分布的一個(gè)重要的考慮因數(shù)。

圖2 未分組數(shù)據(jù)的各種分布與理想可靠度擬合情況Fig.2 Fitting of various distribution and ideal reliability of data ungrouping

4) 將相同數(shù)據(jù)用不同的分布或不同的算法得到擬合曲線與相應(yīng)的可靠度的決定系數(shù)的大小作為擬合好壞的標(biāo)準(zhǔn)還是可行的。如果用3參數(shù)威布爾分布擬合比高斯分布擬合更好的話，那么高斯分布作為一階近似是有道理的。這也是為什么即使是疲勞壽命的數(shù)據(jù)有時(shí)候用高斯分布來擬合也會得到似乎不錯的結(jié)果的原因。

4 威布爾分布的置信區(qū)間：方法1

第3節(jié)討論了用什么標(biāo)準(zhǔn)來判別樣本數(shù)據(jù)用什么分布比較合適或用什么算法得到的參數(shù)更加合適的問題。事實(shí)上數(shù)據(jù)是客觀的，如何去描述是存在主觀意圖的，即從某個(gè)角度去看這些數(shù)據(jù)，以便人們理解和利用這些數(shù)據(jù)。以往各種統(tǒng)計(jì)檢驗(yàn)大都是以高斯分布為基礎(chǔ)，如比較著名的t檢驗(yàn)盡管可以不知道總體的標(biāo)準(zhǔn)差，但仍然是建立在總體分布是高斯分布的基礎(chǔ)之上。所謂置信區(qū)間也是如此，更加重要的是置信區(qū)間僅僅是針對高斯分布的2個(gè)參數(shù)而言的，而且往往是需要已知其中一個(gè)求出另外一個(gè)。問題是對于威布爾分布來說如何求它的置信區(qū)間似乎成了一個(gè)新問題，特別是需要知道各種不同可靠度之下“壽命估計(jì)”，而這種壽命估計(jì)也是需要一個(gè)“置信區(qū)間”的。

那么如何求威布爾分布的置信區(qū)間呢？自然先想到的是求威布爾分布參數(shù)的置信區(qū)間。但同時(shí)求出3個(gè)參數(shù)的置信區(qū)間是相當(dāng)困難的，通過下面的論述將理解這一點(diǎn)，為此要從最小二乘法說起。事實(shí)上，從文獻(xiàn)[1]中推出高鎮(zhèn)同法時(shí)指出3參數(shù)威布爾分布的可靠度為

=1-()=exp{-[(-)]}

(11)

式中:表示可靠度為的疲勞壽命。于是式(11)在取2次自然對數(shù)后可得

ln(ln(1))=ln()-ln()

(12)

式中：=1-/(+1)

再設(shè)

=ln(ln(1/)),=ln(-)

(13)

=+

(14)

式中：

=-ln()，=exp(-)

(15)

故可由計(jì)算機(jī)直接通過最小二乘法求出使得與的相關(guān)系數(shù)最大的，以及相應(yīng)的和的估計(jì)值和，與此同時(shí)也就求出的估計(jì)值了，這就是高鎮(zhèn)同法?，F(xiàn)在的問題是數(shù)據(jù)是在什么程度上服從威布爾分布，盡管通過相關(guān)系數(shù)可認(rèn)為該數(shù)據(jù)與高斯分布相比威布爾分布更合適。但并沒有給出在一定的置信度下該數(shù)據(jù)的分布是否符合威布爾分布的統(tǒng)計(jì)檢驗(yàn)。

注意到此時(shí)的回歸線性方程式(14)的2個(gè)待定系數(shù)與(就是威布爾分布的形狀系數(shù)，而通過及也馬上可得威布爾分布的尺度系數(shù))是通過回歸誤差平方最小化來確定的，服從(0,)分布

=-(+)

(16)

為此可先討論關(guān)于、及的估計(jì)問題。首先要證明由最小二乘法得到的和是最優(yōu)線性無偏估計(jì)量。按最小二乘法

=

(17)

=-

(18)

式中：

(19)

(20)

式中：=(-)，同時(shí)注意到

(21)

這就證明了的線性性。不難證明

(21a)

又因式(18)容易證明也是的線性函數(shù),則

(22)

式中：

=(1-)

(22a)

再注意到式(16)和式(20)，

(23)

在此要注意利用了式(21)和式(21a)。

又因式 (16) 和式(22)可得

(24)

在此要注意到式(22a)

(24a)

由此來證明參數(shù)和的無偏性，即要證明

()=,()=

(25)

事實(shí)上，由式(23)可得

同理可證明()=。

現(xiàn)進(jìn)一步來求的方差

(26)

而

(26a)

在此要注意

1+=()

(26b)

不過至此問題還沒有完全解決，因?yàn)檎`差方差其實(shí)是未知的，因此需要通過可觀察數(shù)據(jù)對其進(jìn)行估計(jì)。注意到可觀察的誤差為

(+)=+(-)+-

(27)

注意到，和的均值都為零

(28)

即可得

(29)

按照假定，

(30)

注意到與的相互獨(dú)立性可得

(30a)

(^)=(-1)+-2=(-2)

(30b)

的無偏估計(jì)為

(31)

這與文獻(xiàn)[4,10]得到的結(jié)果一致。

于是威布爾分布的在置信度=1-(為顯著度)的條件下是否顯著以及威布爾參數(shù)本身的估計(jì)了。從形狀參數(shù)開始。注意到的估計(jì)值是服從高斯分布(,), 其中可由無偏估計(jì)的=/(-2)來代之。

由此可得

=(-)～(-2)

(32)

1) 零假設(shè)，=0(意味著在顯著性水平之下威布爾分布是不正確的)。

2) 備擇假設(shè)，≠0(意味著在顯著性水平之下威布爾分布是正確的)。

||≤2(-2)

(33)

時(shí)接受零假設(shè)，否則接受備擇假設(shè)。

至于的置信區(qū)間就比較容易確定了。此時(shí)用1-來代表置信水平(置信度)，那么又由式(32)可得的置信區(qū)間為

(34)

同理可得關(guān)于的置信區(qū)間，不過此時(shí)，

(35)

注意到式(15)的=-ln()→=exp(-)。 因此的置信區(qū)間為

(36)

于是λ的置信區(qū)間為

exp{[-+2(-2)^]})

(37)

式中

(38)

下面用具體例子來說明。順便指出在此沒有用文獻(xiàn)[4]中的數(shù)據(jù)，是因?yàn)闊o法驗(yàn)證其所得到的結(jié)果，這確實(shí)是很遺憾的事情。

還是利用文獻(xiàn)[2]第136頁中的數(shù)據(jù):{350,380,400,430,450,470,480,500,520,540,550,570,600,610,630,650,670,730,770,840}×10循環(huán)。均值=557，標(biāo)準(zhǔn)差=132.15，中值=545。

同樣可利用Python代碼得到如下結(jié)果：=276.6，=0.999 22，=2.040，=-11.77，=320.98，=0.05,=107.18>_=2.10，則該分布的總體符合威布爾分布?？傮w形狀參數(shù)的置信區(qū)間為(2.020， 2.060)，回歸參數(shù)的置信區(qū)間為(-11.278 -12.265)，總體尺度參數(shù)的置信區(qū)間為(266.2， 385.6)。

由此可見用此法來確定和及的置信區(qū)間的精確度還是相當(dāng)高的，但無法確定安全壽命(位置參數(shù))的置信區(qū)間。此法確定威布爾分布中形狀參數(shù)的置信區(qū)間是有效的，且由此可檢驗(yàn)樣本是在某種置信度下是否符合威布爾分布。雖然可按照線性回歸來給出回歸直線的置信區(qū)間和預(yù)測區(qū)間，但只能是在對數(shù)坐標(biāo)系中進(jìn)行，效果顯然沒有在非對數(shù)坐標(biāo)系中好，第5節(jié)就來解決在某個(gè)可靠度之下疲勞壽命的置信區(qū)間。

5 威布爾分布的置信區(qū)間:方法2

為解決第4節(jié)未解決的問題，需要引入秩分布。事實(shí)上，可設(shè),,…,為相互獨(dú)立且來自同一總體的隨機(jī)變量，且有<<…<, 設(shè)為小于第個(gè)發(fā)生的頻率值，那么容易得到第個(gè)隨機(jī)變量的概率密度函數(shù)為

(39)

即可得

()～Be(,-+1)=

-1(1-)-/(,-+1)

(40)

式中：Be(,)為貝塔分布的概率密度函數(shù)

Be(,)=-1(1-)-1(,)

(40a)

其中：(,)=(+)(()())為貝塔函數(shù)，(·)為伽馬函數(shù)。

為得到置信度為的置信區(qū)間，按文獻(xiàn)[6] 分別設(shè)U,L為置信上限和下限：

(41)

由此的置信區(qū)間為

(42)

式中：≥50%。再由式(41)得

(43)

做變量替換=-+1=1-，注意到的定義式(39)及貝塔分布的定義式(40a)，即式(43)可變?yōu)?/p>

和式(41)的積分限做比較即得

L=1-U(-+1)

(44)

有了式(44)可節(jié)省一半計(jì)算量。在文獻(xiàn)[11]中分別給出了=0.90和=0.95不同的之的表，可以想象在差不多30年前的條件下要計(jì)算這些數(shù)值是非常不容易的事情，最大問題還在于使用是很不方便的，這大概也是這種方法推廣不易的原因之一。而現(xiàn)在計(jì)算機(jī)硬件、軟件如此發(fā)達(dá)，計(jì)算和使用這些數(shù)值完全不需要查表，能夠非常智能化地做到“按需供給”，非常方便人們使用，當(dāng)然算法的基礎(chǔ)還是高鎮(zhèn)同法。

用例3的數(shù)據(jù)，用Python代碼可得，在置信度=0.95、=20的條件下的置信上限:13.91，21.61，28.26，34.37，40.10，45.56，50.78， 55.80，60.64，65.31，69.80，74.13，78.29，82.27，86.04，89.59，92.86，95.78，98.19，99.74。在置信度=0.95、=20 的條件下的置信下限:0.26， 1.81，4.22，7.14，10.41，13.96，17.73，21.71，25.87， 30.20，34.69，39.36，44.20，49.22，54.44，59.90，65.63，71.74，78.39，86.09。按平均秩計(jì)算得到的參數(shù)：=2.040，=320.98，=276.60，=0.999 22， 按上限秩計(jì)算得到的參數(shù)：=2.797，=368.96，=158.80，=0.998 97，按下限秩計(jì)算得到的參數(shù)：=1.986，=343.26，= 333.06，=0.999 26，總體形狀參數(shù)=2.040 的置信區(qū)間為[1.986， 2.797]，總體安全壽命參數(shù)=276.600的置信區(qū)間為[158.800， 333.060]，可靠度=0.999對應(yīng)的安全疲勞壽命= 287.459 置信區(qū)間為[190.019， 343.658]。

在此得到的結(jié)果和文獻(xiàn)[13]中的結(jié)果幾乎是一致的，改進(jìn)了他們的工作，將其智能化了，并且將疲勞壽命的置信曲線都畫出來了，用起來非常方便，且不易搞錯。只要將=0.95換為0.9馬上可得如下結(jié)果：

在置信度=0.9、=20 的條件下置信上限:10.87，18.10，24.48，30.42，36.07，41.49，46.73，51.80，56.73，61.52，66.18，70.71，75.09，79.33，83.41，87.31，90.98，94.36，97.31，99.47。在置信度=0.9、=20的條件下置信下限:0.53，2.69，5.64，9.02，12.69，16.59，20.67，24.91，29.29，33.82，38.48，43.27，48.20，53.27，58.51，63.93，69.58，75.52，81.90，89.13。按平均秩計(jì)算得到參數(shù)：=2.040，=320.98，=276.60，=0.999 22，按上限秩計(jì)算得到參數(shù)：=2.513，=336.90，=204.20，=0.999 05，按下限秩計(jì)算得到參數(shù)：=1.946，=329.56，=327.74，=0.999 31，總體形狀參數(shù)=2.040 的置信區(qū)間為[1.946， 2.513]，總體安全壽命參數(shù)=276.600的置信區(qū)間為[204.200， 327.740]，可靠度=0.999對應(yīng)的安全疲勞壽命=287.459的置信區(qū)間為[225.777， 337.219]。

與=0.95的結(jié)果比較可發(fā)現(xiàn)置信度減少了，那么置信區(qū)間會“縮短”，看起來目標(biāo)收窄了，但命中率(置信度)也小了，可謂“魚與熊掌不可兼得”也，必須權(quán)衡利弊，這也是在一般情況下人們喜歡選取置信度為0.95的原因。圖3和圖4分別為置信度為0.95和0.9的威布爾分布置信限曲線。

圖3 置信度為0.95的疲勞壽命置信限曲線Fig.3 Graph of confidence limit for fatigue life with 0.95 confidence

圖4 置信度為0.9的疲勞壽命置信限曲線Fig.4 Graph of confidence limit for fatigue life with 0.9 confidence

在此有幾個(gè)問題需要強(qiáng)調(diào)一下：

1) 前面所說的置信限是對于分布函數(shù)而言的，而對于可靠度=1-。因此對來講的置信限正好與可靠度的置信限相反，即由的置信上限得到了的置信下限。

2) 對于安全壽命參數(shù)的置信區(qū)間也是與可靠度一致的，即需要“顛倒”。但形狀參數(shù)的置信區(qū)間卻因?yàn)榕c位置參數(shù)是負(fù)相關(guān)的，因此與可靠度相反，反而與的置信限是“一致”的。不難看到在同樣置信水平的條件下由第4節(jié)中得到的關(guān)于的置信區(qū)間要比此節(jié)得到的長度要短，即要好一些。其原因大概是前者是直接計(jì)算得到的，而后者是間接得到的。

3) 另外一個(gè)問題要注意，這里沒有給出尺度參數(shù)的置信區(qū)間，這是因?yàn)橥紶柗植嫉?個(gè)參數(shù)的確定往往是有關(guān)聯(lián)的。先是由相關(guān)系數(shù)極大化來確定出位置參數(shù)，然后再由最小二乘法確定出形狀參數(shù), 最后才間接求出尺度參數(shù)λ。因此得到的所謂其置信區(qū)間就沒有什么意義了(至少用這個(gè)方法是得不到的)。從計(jì)算的結(jié)果也可看到這一點(diǎn)，由平均秩計(jì)算出來的=320.98，居然不在由秩的上下限計(jì)算出來的的區(qū)間[329.56,336.9]里。其實(shí)對于實(shí)際應(yīng)用來說最重要的是由確定的可靠度來定出疲勞壽命的置信區(qū)間。

最后來回顧一下人們對威布爾分布的置信區(qū)間的認(rèn)識是有必要的。由于3參數(shù)數(shù)學(xué)形式的復(fù)雜性，對于3個(gè)參數(shù)的點(diǎn)估計(jì)就存在一定困難，而對于參數(shù)的置信區(qū)間的估計(jì)就難上加難，因此從雙參數(shù)的威布爾分布入手是不奇怪的。但其方法還是比較麻煩，使用起來并不方便，更加重要的是雙參數(shù)忽略威布爾分布中的安全壽命這個(gè)參數(shù)，無論是理論還是實(shí)際應(yīng)用都是欠缺的，當(dāng)然不是沒有人去討論3參數(shù)的置信區(qū)間問題，應(yīng)該說文獻(xiàn)[6]的觀點(diǎn)是正確的，只是使用起來不方便，這一節(jié)利用高鎮(zhèn)同法將其計(jì)算方法加以改進(jìn)，使用起來就沒有問題。文獻(xiàn)[18]的觀點(diǎn)是正確的，但要回到“圖解法”是不可取的。而文獻(xiàn)[19]中給出的3參數(shù)置信區(qū)間所謂近似方法也是相當(dāng)麻煩的，且不說在實(shí)際應(yīng)用中好不好用，最大問題還在于給出的例子是無法“重復(fù)”的。至于文獻(xiàn)[20]例子中的威布爾分布3參數(shù)若用高鎮(zhèn)同法是100%和真值相同，而用該論文的方法還是存在誤差，盡管誤差不算大，其他例子也存在同樣問題。其提出用自助法來優(yōu)化置信區(qū)間的思想是可取的，但沒有講清楚如何得到得到威布爾3參數(shù)點(diǎn)估計(jì)的具體方法，而這恰恰是最基本的。

6 結(jié)束語

將大樣本數(shù)據(jù)進(jìn)行分組實(shí)在是某種“無奈之舉”，其中一個(gè)目的在于檢驗(yàn)數(shù)據(jù)是否服從正態(tài)分布，過程相當(dāng)麻煩，而且一旦檢驗(yàn)出不服從，就沒有下文了。在疲勞統(tǒng)計(jì)學(xué)中研究的數(shù)據(jù)若不服從正態(tài)分布，那么在一般情況下往往是符合威布爾分布。此時(shí)可用高鎮(zhèn)同法通過相關(guān)系數(shù)或決定系數(shù)的大小來判定樣本數(shù)據(jù)服從威布爾分布是否更合適，且根本無需將數(shù)據(jù)“分組”。

當(dāng)然不服從威布爾分布也可能是屬于其他類型的分布，可根據(jù)第4節(jié)中給出的方法做一個(gè)關(guān)于形狀參數(shù)的t檢驗(yàn)，以確定在一定的置信度該數(shù)據(jù)是否屬于威布爾分布。

由于威布爾分布的3個(gè)參數(shù)具有內(nèi)在的某種復(fù)雜的非線性關(guān)系，因此很難同時(shí)給出3個(gè)參數(shù)的置信區(qū)間(這個(gè)問題將在另一篇論文中解決)。而只能提供可靠度的上下限置信曲線，實(shí)際上是提供了在一定可靠度下疲勞壽命的置信區(qū)間，而這個(gè)置信區(qū)間恰恰是人們在解決可靠性問題中最想知道的信息！這也是一個(gè)智能化解決方案，其基礎(chǔ)仍然是高鎮(zhèn)同法。

從更深層次的角度看，一組數(shù)據(jù)或者說一個(gè)樣本，從本質(zhì)上講是該數(shù)據(jù)的總體在某個(gè)參照系(維度)上的投影。利用各種統(tǒng)計(jì)方法就是希望從這個(gè)投影中找到其總體的廬山真面目。如果總體比較簡單，比如就是一個(gè)簡單的高斯分布，投影又沒有受到“扭曲”，那么就比較容易判別出其總體的高斯分布的2個(gè)主要參數(shù)。但如果不是高斯分布，或者即使是高斯分布投影受到扭曲，比如數(shù)據(jù)不夠多，或者只是投影“片面”的數(shù)據(jù)，那么就很難得出關(guān)于其總體屬于何種分布的結(jié)論。其實(shí)所謂高斯分布也好、威布爾分布也好或者其他分布也好，都是人類在經(jīng)驗(yàn)中總結(jié)出來的分布，是人們用來描述世界的工具。在許多場合中是有效的，但極有可能在某些場合中是不那么有效甚至無效的，即不存在無條件是正確的東西。統(tǒng)計(jì)的任務(wù)是在不確定中尋找一定條件下的某種確定性，某種能夠量化的東西，那就是“概率”。本文也是一樣，無法提供區(qū)分高斯分布和威布爾分布的絕對標(biāo)準(zhǔn)，只能是提供與理想可靠度決定系數(shù)的量化標(biāo)準(zhǔn)，在給定的置信度和可靠度之下疲勞壽命的置信區(qū)間。

文獻(xiàn)[1]和本文取得一些成果僅僅是疲勞統(tǒng)計(jì)學(xué)智能化的一個(gè)“開始”，還處于“初級階段”。希望能起到拋磚引玉的作用，在不遠(yuǎn)的將來能看到取得更多的智能化成果。

致謝

感謝萬偉浩先生對有關(guān)研究工作的全力支持。感謝鮑蕊教授無私的幫助和指導(dǎo)。