基于EM算法的Marshall-Olkin二元指數分布的參數估計

袁守成，張 賓，陳相兵

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

0 引言

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

1 預備知識

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

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

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

因而，X1和X2的聯合概率密度為:

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

2 參數估計的EM算法

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

對式(1)分別關于 λ1，λ2，λ3求導，有：

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

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

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

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

由于:

因此:

類似地，計算可得:

從而，偽似然函數可以寫為:

上式簡寫為：

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

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

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

3 區間估計

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

其中，zα為標準正態分布的上側α分位點。

4 數值模擬

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

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

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