基于EM算法的Dirichlet分布的參數估計

吳 琴，葛文秀

(華南師范大學數學科學學院，廣東 廣州 510631)

0 引言

成分數據是記錄觀測數據相對關系的多維數據，在科學研究和日常生活中經常用到，比如地質學中巖石的組成成分比例、家庭消費模式中不同類別所占百分比等.Dirichlet分布是處理成分數據的重要工具，但是由于隨機變量各分量之間非負定和的約束限制，其統計分析具有較大的難度，傳統的參數估計方法算法復雜且不能保證有效性，極大地限制了其統計分析的進展.許多學者致力于Dirichlet分布的研究.張堯庭[1]詳細給出了Dirichlet分布[2]參數估計的算法，但是其提出的算法由于計算復雜且不能保證是有效估計(非負參數的估計值為負值)從而限制了其應用.Hijazi等[3]給出了如何計算Dirichlet回歸模型參數的極大似然估計，但是其使用的Newton-Raphson算法嚴重依賴初值的選取而且存在不收斂的缺陷.Ma[4]提出了Dirichlet分布的Bayes估計，但有時存在參數不可估的現象.為了克服上述缺點，得到Dirichlet分布參數的有效估計，本文將根據Dirichlet分布的性質提出其隨機表示[5]，通過隨機表示的具體形式引入缺失數據，構造Dirichlet分布參數估計的EM算法[6].本文提出的EM算法，不僅不依賴于初值的選取，而且能保證得到的估計為正數，也就是有效估計.最后，統計模擬結果表明，本文方法得到的估計精度良好，即使樣本量不大(比如n=20)，估計的偏差也是可以接受的，當樣本量n=500的時候，估計值非常接近真實值.

1 預備知識

定義1[5]稱m維隨機向量X=(X1，…，Xm)T服從參數是α=(α1，…，αm)T的Dirichlet分布，記為X～Dirichlet(α)，如果其概率密度函數如下：

其對數似然函數為

2 理論和方法

2.1 理論基礎

從引理1中可以看出，Gamma分布中的參數β在Dirichlet分布中消失了.基于引理1，可以構建如下隨機表示：

證明令引理1中的β=1即可.

引理2的證明可參考文獻[1].

利用引理2的結論做變換Yj=sxj即可得到引理3的結論，變換的雅可比行列式為1/xj.

定理1 若X=(X1，…，Xm)T服從參數為α=(α1，…，αm)T的Dirichlet分布，X為觀測數據，Y=(Y1，…，Ym)T為X的基向量，即Yj～Gamma(αj，1)，j=1，…，m.則

證明

2.2 Dirichlet參數的極大似然估計

EM算法的M步為求解完全數據的對數極大似然關于參數求偏導的方程組：

但上述方程組沒有顯式解，用Newton-Raphson迭代法求解，具體迭代方法為

E步為求基于觀測數據下缺失數據logyij的條件期望，由定理1得：

3 統計模擬

表1 不同參數下EM算法估計的模擬結果

從表1的結果中可以看出，即使樣本量很小(n=20)，EM算法得到的估計的偏差也是可以接受的，當樣本量n=100的時候，估計的誤差已經很小了，當n=500的時候，估計值非常接近真實值.此外，對于m=2和m=3兩種情況對比，發現估計的功效差別不大.實際上，本文嘗試了m取其他值的情況，比如m=5，其模擬結果顯示功效與m=2時并無明顯差異.可見本文所提出的EM方法估計精度不會隨著m的增大而降低.