吳 琴,葛文秀
(華南師范大學數學科學學院,廣東 廣州 510631)
成分數據是記錄觀測數據相對關系的多維數據,在科學研究和日常生活中經常用到,比如地質學中巖石的組成成分比例、家庭消費模式中不同類別所占百分比等.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[5]稱m維隨機向量X=(X1,…,Xm)T服從參數是α=(α1,…,αm)T的Dirichlet分布,記為X~Dirichlet(α),如果其概率密度函數如下:
其對數似然函數為
從引理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.則
證明
EM算法的M步為求解完全數據的對數極大似然關于參數求偏導的方程組:
但上述方程組沒有顯式解,用Newton-Raphson迭代法求解,具體迭代方法為
E步為求基于觀測數據下缺失數據logyij的條件期望,由定理1得:
表1 不同參數下EM算法估計的模擬結果
從表1的結果中可以看出,即使樣本量很小(n=20),EM算法得到的估計的偏差也是可以接受的,當樣本量n=100的時候,估計的誤差已經很小了,當n=500的時候,估計值非常接近真實值.此外,對于m=2和m=3兩種情況對比,發(fā)現估計的功效差別不大.實際上,本文嘗試了m取其他值的情況,比如m=5,其模擬結果顯示功效與m=2時并無明顯差異.可見本文所提出的EM方法估計精度不會隨著m的增大而降低.