李丹青， 熊文潔， 張正成

(1.蘭州交通大學(xué) 數(shù)理學(xué)院，甘肅 蘭州 730070；2.海南師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，海南 ?？?571158)

在可靠性理論中，Rosen[1]率先使用負(fù)載共享系統(tǒng)模型來(lái)研究復(fù)合材料的可靠性，介紹了一種復(fù)合材料失效的理論和實(shí)驗(yàn)方法，一束纖維可視為受恒定負(fù)載影響的并聯(lián)系統(tǒng)，當(dāng)受到平行于纖維方向的單軸拉伸負(fù)載時(shí)，纖維被視為一個(gè)統(tǒng)計(jì)分布，在施加的應(yīng)力下導(dǎo)致纖維失效。之后，關(guān)于可靠性系統(tǒng)的研究便滲入了各行各業(yè)，統(tǒng)計(jì)學(xué)家們紛紛開(kāi)始研究可靠性系統(tǒng)并將其運(yùn)用到實(shí)際中。Boland 等人[2]討論了在系統(tǒng)中如何分配冗余組件的位置從而優(yōu)化系統(tǒng)壽命的問(wèn)題，他們考慮了主動(dòng)冗余和備用冗余，特別是在串聯(lián)和并聯(lián)系統(tǒng)中獲得了廣泛的結(jié)果。在大多數(shù)關(guān)于負(fù)載共享系統(tǒng)的研究中，系統(tǒng)的性能僅針對(duì)元件的壽命遵循指數(shù)分布的情況進(jìn)行了研究，Yun 和Cha[3]研究了當(dāng)系統(tǒng)中元件的壽命是任何連續(xù)的隨機(jī)變量時(shí)，考慮負(fù)載共享并聯(lián)系統(tǒng)的可靠性及其負(fù)載分配問(wèn)題。關(guān)于負(fù)載共享系統(tǒng)的可靠性更詳細(xì)的問(wèn)題可見(jiàn)文獻(xiàn)[4-9]。

Kumaraswamy分布[10]是在1980 年由Kumaraswamy研究水力學(xué)當(dāng)中的日常降雨量提出的，它是一種在(0，1) 范圍內(nèi)取值的雙參數(shù)連續(xù)型分布，與Beta 分布有許多相同的特征但其分布函數(shù)有非常簡(jiǎn)單的顯式，這個(gè)特征對(duì)這個(gè)分布的模擬計(jì)算與應(yīng)用帶來(lái)很大的方便。Nadarajah[11]給出了Kumaraswamy分布的幾種一般形式及其相關(guān)應(yīng)用。MITNIK[12-13]在過(guò)去研究的基礎(chǔ)上擴(kuò)展了Kumaraswamy分布的屬性特征，并提出了兩種簡(jiǎn)單的Kumaraswamy分布變量屬性：線(xiàn)性轉(zhuǎn)化和取冪形式下該分布具有緊湊型特征.該分布作為壽命模型在可靠性方面非常有用，但關(guān)于負(fù)載共享并聯(lián)系統(tǒng)的研究卻被忽略了。

Kumaraswamy分布的可靠性函數(shù)：

(1)

Kumaraswamy分布的密度函數(shù)：

f(x)=αβxα-1(1-xα)β-1,00,β>0

(2)

在早期，KIM 和KVAM[14]考慮了k個(gè)元件并聯(lián)的負(fù)載共享系統(tǒng)，元件的故障率取決于系統(tǒng)中剩余存活的元件個(gè)數(shù).起初，各元件獨(dú)立同分布且具有共同的失效率θ，直到第一個(gè)元件失效后，剩余存活元件的失效率為γ1θ，第二個(gè)元件失效后，剩余存活元件的失效率為γ2θ(γi>0,i=1,2,…,k-1)等等。在基于負(fù)載共享規(guī)則是未知的前提下，根據(jù)最大似然估計(jì)導(dǎo)出負(fù)載共享系統(tǒng)參數(shù)γi,θ的統(tǒng)計(jì)推斷方法。KVAM 和PENA[15]基于獨(dú)立同分布的k個(gè)元件并聯(lián)的負(fù)載共享系統(tǒng)，導(dǎo)出了等負(fù)載分配模型中負(fù)載分配參數(shù)的估計(jì)值。在等負(fù)荷分配模型中，在第一個(gè)元件失效后，其余元件的故障率從r(t)變?yōu)棣?r(t)，然后在下一次故障后變?yōu)棣?r(t)等等.然后，他們獲得了元件基線(xiàn)累積危險(xiǎn)函數(shù)的半?yún)?shù)估計(jì)，并討論了它的漸近分布。SINGH 等[16]通過(guò)在開(kāi)始時(shí)使用恒定負(fù)載，然后在一定數(shù)量的元件故障后線(xiàn)性增加失效率，進(jìn)一步修改了KIM 和KVAM[14]的模型。PARK[17]研究了當(dāng)元件的壽命分布是指數(shù)或Weibull 分布時(shí)這類(lèi)系統(tǒng)的參數(shù)，并且在元件為指數(shù)分布的情況下，采用KIM 和KVAM[14]中的模擬數(shù)據(jù)，在相等負(fù)載共享的規(guī)則下獲得了閉合形式的MLE 和漸近最佳無(wú)偏估計(jì)(ABUE) .HUANG[18]研究了EKD 家族(Exponentiated Kumaraswamy-Dagum) 分布函數(shù)性質(zhì)，包括密度函數(shù)、風(fēng)險(xiǎn)函數(shù)、可靠性等，并引入實(shí)例比較EKD 與其他函數(shù)的差異.最近，ZHANG 和Balakrishnan[19]研究了一種通用的負(fù)載共享并聯(lián)系統(tǒng)，其中系統(tǒng)元件的壽命是任意連續(xù)的隨機(jī)變量，然后考慮具有兩個(gè)獨(dú)立元件的串聯(lián)系統(tǒng)中一個(gè)負(fù)載備用的最佳分配問(wèn)題，并討論了某些特定系統(tǒng)的參數(shù)的最大似然估計(jì)。

1 建立模型

1.1 模型的描述

首先我們建立一個(gè)負(fù)載共享并聯(lián)系統(tǒng)的模型來(lái)估計(jì)出它的參數(shù)，假設(shè)這個(gè)模型具有以下特征:

⑴n個(gè)獨(dú)立同分布的負(fù)載共享系統(tǒng)中各包含k個(gè)元件，它們帶有一個(gè)共同的恒定負(fù)載,該負(fù)載在其工作元件中均勻分配;

⑵對(duì)這n個(gè)系統(tǒng)進(jìn)行測(cè)試并記錄其故障時(shí)間，設(shè)Tij(i=1,2,…,n;j=1,2,…,k)是第i個(gè)并聯(lián)系統(tǒng)中第(j-1)和第j個(gè)元件失效之間的時(shí)間間隔;

⑶系統(tǒng)內(nèi)一個(gè)元件的失效會(huì)增加剩余存活元件的負(fù)載;

⑷起初，元件的壽命服從參數(shù)為α1,β1的Kumaraswamy分布，在元件第一次失效后，初始參數(shù)由α1,β1變?yōu)棣?,β2，第二次失效后，參數(shù)由α2,β2變?yōu)棣?,β3，直到元件第(k-1)次失效后，最后一個(gè)元件的參數(shù)變?yōu)棣羕,βk;

⑸負(fù)載共享參數(shù)α1,α2,…,αk和β1,β2,…,βk是獨(dú)立且未知的。

1.2 系統(tǒng)密度函數(shù)的建立

從上述系統(tǒng)模型的特征中可以看出，起初，在具有k個(gè)元件的系統(tǒng)中，其失效率為h(t,α1,β1);當(dāng)元件第一次失效后，系統(tǒng)剩余存活的(k-1)個(gè)元件的失效率為h(t,α2,β2);在第二次失效后，系統(tǒng)剩余存活的(k-2)個(gè)元件的失效率為h(t,α3,β3);最后，當(dāng)?shù)?k-1)個(gè)元件失效后，最后存活的那個(gè)元件的失效率為h(t,αk,βk)。

其中，服從Kumaraswamy分布的單個(gè)元件的失效率函數(shù)為：

h(t,α,β)=αβxα-1(1-xα)-1,00,β>0

(3)

那么，在(j-1)個(gè)元件完全失效后第j個(gè)元件的條件失效率函數(shù)為：

hTj|T(j-1)(t)=(k-j+1)h(t,αj,βj),j=1,2,…,k

(4)

根據(jù)以上假設(shè)(1)-(5) 及公式(4) 可以得出，隨機(jī)變量Tij服從具有失效率(k-j+1)×h(t,αj,βj)的Kumaraswamy分布。所以，在第i個(gè)系統(tǒng)中第j個(gè)元件失效時(shí)間的密度函數(shù)為：

f(tij)=(k-j+1)αjβjtijαj-1·(1-tijαj)(k-j+1)βj-1，00,β>0

(5)

2 參數(shù)估計(jì)

設(shè)A=(α1,α2,…,αk)′，B=(β1,β2,…,βk)′存在于空間(0,∞)k中，其中αj>0，βj>0(j=1,2,…,k)，則第i個(gè)負(fù)載共享并聯(lián)系統(tǒng)模型的似然函數(shù)為：

基于n個(gè)獨(dú)立同分布系統(tǒng)的似然函數(shù)為：

這里，T={tij,i=1,2,…,n;j=1,2,…,k}。所以，相應(yīng)的對(duì)數(shù)似然函數(shù)為：

然后便可以根據(jù)這相應(yīng)的2k個(gè)對(duì)數(shù)似然方程：

(6)

(7)

(8)

(9)

像高斯-賽德?tīng)柕ň湍芎芎玫厍蠼獬鲞@2k個(gè)參數(shù)的值,一般分為以下幾步：

⑴選擇合適的初始值α1(0),…,αk(0)，將其代入Nj=0中求出βk的解，將其表示為βk(1)；

⑵將βk(1)代入Mj=0中求出αk的解，將其表示為αk-1(1)；

⑶同樣將α1(0),…,αk-1(0),αk(1)代入Nj=0中求出βk-1的解，將其表示為βk-1(1)；

⑷將βk-1(1)代入Mj=0中求出αk的解，將其表示為αk-1(1)；

⑸依次類(lèi)推，在最后的表達(dá)式中固定其他變量來(lái)求解αj和βj，得出αj(1),βk-1(1)使得Mj=0,Nj=0；

接下來(lái)，為簡(jiǎn)單起見(jiàn)，我們來(lái)研究α=1時(shí)Kumaraswamy分布的最大似然估計(jì)值(MLE) 。此時(shí)，具有一個(gè)參數(shù)β的可靠性函數(shù)是：

相應(yīng)的概率密度函數(shù)為：

f(x)=β(1-x)β-1,0x1,β>0

并且第i個(gè)系統(tǒng)的第j個(gè)失效時(shí)間的概率密度函數(shù)為：

f(tij|αj)=(k-j+1)βj(1-tij)(k-j+1)βj-1，i=1,2,…,n;j=1,2,…,k。

(10)

所以，第i個(gè)系統(tǒng)(i=1,2,…,n)的似然函數(shù)是：

且基于n 個(gè)樣本的似然函數(shù)是：

(11)

則其對(duì)數(shù)似然函數(shù)為：

通過(guò)對(duì)上式進(jìn)行微分，我們可以得出一個(gè)對(duì)數(shù)似然方程：

因此，我們便可以得出其最大似然估計(jì)值(MLE) 為：

根據(jù)相對(duì)速率(RE) ，我們可以得出ABUE 優(yōu)于MLE ，即使在n→∞時(shí)它的極限趨于1 ，但是這細(xì)微的差異對(duì)于小樣本來(lái)說(shuō)也是不可忽視的。此外，

顯然，k×kHessian 矩陣{?2l/?βi?βj}是負(fù)定的，所以在這個(gè)子空間中存在一個(gè)向量B使 得產(chǎn)生(11) 的全局最大值。

在空間(0,∞)k中讓?duì)?(β1,β2,…,βk)′，這里，βj>0(j=1,2,…,k)。那么，參數(shù)向量Β的Fisher 信息矩陣為

為參數(shù)構(gòu)建基于似然比的聯(lián)合置信度，觀察到的Fisher 信息矩陣F的倒數(shù)提供了大樣本正態(tài)分布中的協(xié)方差的估計(jì)。對(duì)于大樣本，Β∈(0,∞)k的近似1-α置信橢球?yàn)?/p>

3 模擬研究

圖1 (β1,β2)的(80%，90%，95%，99%) 的置信橢圓圖

表1 失效時(shí)間數(shù)據(jù)

表2 MLE 和ABUE

表3 MLE和ABUE的經(jīng)驗(yàn)偏差(bias)及均方誤差(MSE)

在下文中，分別構(gòu)造了基于Bootstrap 技術(shù)的兩種參數(shù)置信區(qū)間，如百分位bootstrap(boot-p) 和bootstrap-t(Boot-t) 方法，更多細(xì)節(jié)詳見(jiàn)Efron[21]和Hall[22]。

圖2 (β1,β3)的(80%，90%，95%，99%) 的置信橢圓圖

3.1 Boot-p算法

⑴用公式(10) 生成n個(gè)樣本Tij(i=1,2,…,n;j=1,2,…,k)；

⑶重復(fù)步驟2B次；

3.2 Boot-t 算法

⑴用公式(10) 生成n個(gè)樣本Tij(i=1,2,…,n;j=1,2,…,k)；

⑷重復(fù)步驟2-3B次；

在表4中，我們?nèi)耘f設(shè)置k=3,β1=1,β2=2,β3=5。并且在R 語(yǔ)言軟件的幫助下，我們可以在樣本數(shù)n分別為20，50，75，100的情況下計(jì)算不同的估計(jì)值，標(biāo)準(zhǔn)誤差(SE) 及 其95% Bootstrap 置信區(qū)間(CI) 的寬度.這里我們采用了5000次Bootstrap程序(B=5000)。從表4可以看出，當(dāng)樣本數(shù)逐漸增多時(shí)，參數(shù)估計(jì)的標(biāo)準(zhǔn)誤差越小，區(qū)間寬度越小，故參數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)更精確，且由于Boot-t 是由Boot-p 方法進(jìn)行改進(jìn)得到的，因而比Boot-p 方法更加準(zhǔn)確.結(jié)合表2 和表4 我們可以看出，當(dāng)樣本量同為20 時(shí)，Bootstrap 估計(jì)值均優(yōu)于MLE 和ABUE ，由此可見(jiàn)，Bootstrap 技術(shù)的優(yōu)越性。

表4 95% bootstrap 置信區(qū)間(CI)

4 總結(jié)

在本文中，我們研究了負(fù)載共享系統(tǒng)可靠性的參數(shù)估計(jì)方法.在壽命遵循一般的負(fù)載 共享原則的基礎(chǔ)上，我們考慮了元件的基礎(chǔ)分布為kumaraswmy分布。

我們給出了適當(dāng)?shù)牡椒▉?lái)討論具有兩個(gè)未知參數(shù)的分布的參數(shù)估計(jì)器。此外，構(gòu)建了估計(jì)統(tǒng)計(jì)量的漸近分布。對(duì)于簡(jiǎn)化的模型，我們通過(guò)模擬的數(shù)據(jù)提供了閉合形式的MLE，ABUE 和參數(shù)的置信區(qū)間。