袁守成，張 賓，陳相兵

（1.普洱學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，云南 普洱 665000；2.湖北民族學(xué)院 教育學(xué)院，湖北 恩施 445000；3.四川大學(xué)錦江學(xué)院，成都 620860）

0 引言

多元指數(shù)分布在工業(yè)技術(shù)上有重要的應(yīng)用，基于不同的系統(tǒng)可靠性可以導(dǎo)出不同的多元指數(shù)分布形式。在眾多形式中，Marshall-Olkin多元指數(shù)分布[1]是一類唯一保持“無記憶性”特征的分布，其中二元隨機(jī)變量結(jié)構(gòu)更是被廣泛關(guān)注。Marshall-Olkin二元指數(shù)分布是既不被計數(shù)測度控制，又不被Lebesgue測度控制的分布[2]，并且與一元指數(shù)分布不同，Marshall-Olkin二元指數(shù)分布不屬于指數(shù)族，其概率密度函數(shù)具有奇異部分，標(biāo)準(zhǔn)的參數(shù)估計方法已不再適用，需要用新的處理方法和手段對其參數(shù)進(jìn)行研究。Arnold(1968)[3]用矩方法對其參數(shù)做了研究，并得到了多元指數(shù)分布參數(shù)的無偏估計；李國安(1993)[4]從二元指數(shù)分布的識別性[5]入手，利用極大似然估計(MLE)法討論了參數(shù)的估計問題，并獲得該分布參數(shù)的一致最小方差無偏估計；進(jìn)一步，李國安(2015)[6]討論了指數(shù)分布的抽樣定理，并給出了四參數(shù)Marshall-Olkin二元指數(shù)分布參數(shù)的一致最小方差無偏估計。本文基于EM算法，通過引入新的潛在變量，克服了直接利用MLE確定參數(shù)的困難，給出了Marshall-Olkin二元指數(shù)分布的參數(shù)的點估計和區(qū)間估計，通過數(shù)值模擬，發(fā)現(xiàn)該方法切實可行。

1 預(yù)備知識

設(shè)某隨機(jī)變量U服從指數(shù)分布，那么它的密度函數(shù)和生存函數(shù)可分別表示為：

其中 Ι(0，∞)(u)表示 u 在 (0，∞)上示性函數(shù)，λ＞0 為尺度參數(shù)，記作U～E(λ)。假設(shè)三個相互獨立的隨機(jī)變量U1，U2，U3分 別 有 U1～E(λ1) ，U2～E(λ2) ，U3～E(λ3) ，若X1=min(U1，U3)，X2=min(U2，U3)，則二元隨機(jī)向量 (X1，X2)服從參數(shù)為 λ1，λ2，λ3的Marshall-Olkin二元指數(shù)分布，記作(X1，X2)～MOBE(λ1，λ2，λ3)。

若 (X1，X2)～MOBE(λ1，λ2，λ3)，令 z=max{x1，x2} ，則它的生存函數(shù)為：

因而，X1和X2的聯(lián)合概率密度為:

其中，f1(x1，x2)=f(x1;λ1)f(x2;λ2+λ3)，f2(x1，x2)=f(x1;λ1+λ3)

2 參數(shù)估計的EM算法

設(shè)有來自總體 MOBE(λ1，λ2，λ3)樣本容量為 n 的觀測值 (x11，x21)，…，(x1n，x2n)，可以做如下標(biāo)記：I1={i:x1i＜x2i}，I2={i:x1i＞x2i},I3={i:x1i=x2i=x},|I1|=n1,|I2|=n2，|I3|=n3，其中 ||Ii表示集合Ii的元素個數(shù)，那么MOBE的對數(shù)似然函數(shù)可以寫為：

對式(1)分別關(guān)于 λ1，λ2，λ3求導(dǎo)，有：

容易發(fā)現(xiàn),參數(shù) λ1，λ2，λ3的 MLE 無法用顯式的形式給出來,而必須要通過求解非線性方程組才能得到,這給參數(shù)的估計帶來不便。下面基于EM算法對參數(shù) λ1，λ2，λ3的MLE給予改進(jìn)。

考慮到存在著一個與二元隨機(jī)向量(X1，X2)相聯(lián)系的潛在二元隨機(jī)向量 (Δ1，Δ2),其表示為:

可以看到,就算 (X1，X2)已知時,(Δ1，Δ2)也未必知道。若 X1=X2,則有 Δ1=Δ2=0,這是確定的,但是如果X1≠X2,那么 (Δ1，Δ2)就不確定了。例如,當(dāng) (x1，x2)∈I1時,(Δ1，Δ2)的可能取值為 (1，0)或 (1，2)。類似地,當(dāng) (x1，x2)∈I2時,(Δ1，Δ2)的可能取值為 (0，2)或 (1，2)。

為了實現(xiàn)EM算法，本文采用類似于文獻(xiàn)[7,8]中的方法。在“E步”，如果 (X1，X2)∈I3，那么 (Δ1，Δ2)是完全已知的，這時對數(shù)似然函數(shù)(1)的貢獻(xiàn)不變；如果 (X1，X2)∈I1或(X1，X2)∈I2，把它們當(dāng)做缺失數(shù)據(jù)處理。當(dāng) (x1，x2)∈I1時,(Δ1，Δ2)的可能取值為 (1，0)或 (1，2)，認(rèn)為對數(shù)似然函數(shù)的貢獻(xiàn)是由 (x1，x2，u1)和 (x1，x2，u2)兩部分形成，其中 u1和 u2分別為已知 (x1，x2)∈I1時 (Δ1，Δ2)取 (1，0)和 (1，2)的條件概率，而這時的對數(shù)似然函數(shù)也被稱為偽對數(shù)似然函數(shù)。類似 地 ，當(dāng) (x1，x2)∈I2時,(Δ1，Δ2) 的 可 能 取 值 為 (0，2) 或(1，2)，認(rèn)為偽對數(shù)似然函數(shù)的貢獻(xiàn)由 (x1，x2，v1)和 (x1，x2，v2)兩部分形成，其中 v1和 v2分別為已知 (x1，x2)∈I2時 (Δ1，Δ2)取(0，2)和(1，2)的條件概率．在“M步”中，對偽對數(shù)似然函數(shù)關(guān)于未知參數(shù)求導(dǎo)，從而確定最大值。下面計算條件概率 u1，u2，v1，v2。

由于:

因此:

類似地，計算可得:

從而，偽似然函數(shù)可以寫為:

上式簡寫為：

在“M 步”，對式(3)極大化，分別對 λ1，λ2，λ3求導(dǎo)，可以得到:

下面給出EM算法的具體實施步驟：

（1）選 取 初 值 ：給 出 參 數(shù) λ1，λ2，λ3的 初 始 估 計 值

3 區(qū)間估計

為表述方便，用V=(vij)1≤i，j≤3表示 Fisher觀測信息矩陣 I-1(θ0)，則參數(shù) λi，i=1，2，3在置信水平為1-α 下的置信區(qū)間為：

其中，zα為標(biāo)準(zhǔn)正態(tài)分布的上側(cè)α分位點。

4 數(shù)值模擬

為了驗證EM算法對MOBE分布參數(shù)估計的理論結(jié)果，本文通過數(shù)值模擬說明其估計性能。假定U1，U2，U3分布的真實參數(shù)為λ1=λ2=λ3=1，分別取樣本量n＝25,50,100,150,200。由于EM算法是收斂算法，故估計的結(jié)果與初始值的取值大小無關(guān)，不妨設(shè)定初始值，并設(shè)定δ=10-5。用以下指標(biāo)對EM算法的估計性能進(jìn)行考察，分別為平均估計（記為AE），均方誤差（記為MSE）以及置信水平為95%的區(qū)間覆蓋率（記為CP），模擬次數(shù)設(shè)定為1000次。

表1 MOBE分布中參數(shù) λ1，λ2，λ3的AE,MSE及CP

通過表1，容易發(fā)現(xiàn)當(dāng)樣本量小的時候，λ1，λ2，λ3的點估計和區(qū)間估計會有一些小偏差，但是隨著樣本量的增大，估計越來越接近真實值，均方誤差也越來越小，同時區(qū)間覆蓋率也越來越精確。從而，用EM算法給出的估計是一種相合的估計，使用該方法估計MOBE中參數(shù)是一種行之有效的方法。