鄭倩貞， 徐平峰

(長春工業(yè)大學 數(shù)學與統(tǒng)計學院， 吉林 長春 130012)

0 引 言

在高斯圖模型中，協(xié)方差逆陣可解釋隨機變量間的條件獨立性。矩陣中的零元素表示對應兩個隨機變量間條件獨立，因此,協(xié)方差逆陣中零元素的估計問題是圖模型選擇中的關鍵問題。懲罰似然方法可同時實現(xiàn)模型選擇和參數(shù)估計，其中Lasso懲罰是常用的求解稀疏估計的懲罰項，例如：Yuan M等[1]、Banerjee O等[2]、Friedman J等[3]、Stidler N等[4]、徐平峰等[5]都采用l1懲罰對數(shù)似然方法估計協(xié)方差逆陣,從而得到稀疏的圖結構。文獻[5]基于gCoda方法將成分數(shù)據的圖模型選擇問題轉換為高斯圖模型的協(xié)方差逆陣估計問題；文獻[4]研究的是數(shù)據存在隨機缺失時的高斯圖模型選擇問題；而其余三者關注的均是數(shù)據完全觀測下的高斯圖模型選擇問題，但三者采用的優(yōu)化算法不同。

研究表明，Lasso懲罰可得到稀疏估計，但其估計結果會產生偏差[6]。Fan J等[6]提出的SCAD懲罰和Zou H[7]提出的自適應Lasso懲罰都能夠降低估計的偏差。Fan J等[8]將這兩種懲罰應用于完全數(shù)據的圖模型選擇問題以減弱協(xié)方差逆陣估計的偏差。但在實際高維復雜數(shù)據中，往往存在數(shù)據部分缺失的情況。針對數(shù)據完全隨機缺失情形，文中將采用自適應Lasso懲罰對數(shù)似然方法對高斯圖模型選擇進行研究，并與Stidler N等[1]的l1懲罰對數(shù)似然方法(MissGLasso)進行比較。

1 自適應MissGLasso

設x=(xobs,xmis)為來自多元高斯分布Np(μ,Ω-1)的n個獨立同分布樣本。令樣本中的觀測數(shù)據xobs=(xobs,1,…,xobs,n)，缺失數(shù)據xmis=(xmis,1,…,xmis,n)，其中xobs,i和xmis,i分別表示第i個樣本的觀測數(shù)據和缺失數(shù)據的集合，i=1,2,…,n。與文獻[4]相同，文中假設數(shù)據的缺失機制為可忽略的?；谟^測數(shù)據xobs，可給出觀測對數(shù)似然函數(shù)

l(μ,Ω;xobs)=

((Ω-1)obs,i)-1×(xobs,i-μobs,i))，

(1)

式中：μobs,i----第i個樣本中觀測變量的均值;

(Ω-1)obs,i----第i個樣本中觀測變量的協(xié)方差陣。

在懲罰似然的框架下，下述優(yōu)化問題的解可作為均值μ和協(xié)方差逆陣Ω的估計，

(2)

式中：λ----正則化參數(shù);

ωjk----Ω中的(j,k)元素;

p(·)----定義在各元素ωjk上的懲罰函數(shù)，j,k=1,2,…,p。

若令懲罰函數(shù)p(|ωjk|)=|ωjk|，則可得到Lasso懲罰，那么由上式求得的解即為MissGLasso估計[4]。

(3)

我們稱該估計為自適應MissGLasso(AdaMissGLasso)。針對上述優(yōu)化問題，可結合EM算法[9]和GLasso算法[3]進行優(yōu)化求解。

2 EM算法

EM算法常用于處理缺失數(shù)據下的參數(shù)估計問題，該算法通過迭代交替進行E(期望)步和M(最大化)步，直至算法收斂。E步是在給定觀測數(shù)據和當前參數(shù)的情況下，計算完全數(shù)據的對數(shù)似然函數(shù)的條件期望；M步是關于待估參數(shù)極大化E步中的完全似然條件期望。令第t次迭代的估計值(μ(t),Ω(t))作為當前參數(shù)，其中(μ(0),Ω(0))代表EM算法的初始值，設為均值插補后的GLasso估計。第(t+1)次迭代如下:

1)E步。

給定觀測數(shù)據xobs和當前參數(shù)(μ(t),Ω(t))，計算完全數(shù)據懲罰對數(shù)似然函數(shù)的條件期望，記為Q函數(shù)。為簡單起見，記E(·|xobs,μ(t),Ω(t))為Et(·)。于是，Q函數(shù)可表示為

Q(μ,Ω|μ(t),Ω(t))=

(4)

其中

與Fan J等[8]相同，文中取γ=0.5。

(5)

因此，E步只需計算Et(xij)和Et(xijxik)。

對第i個樣本xi=(xobs,i,xmis,i)，μ和Ω可分塊表示為:

給定觀測數(shù)據xobs,i和當前參數(shù)(μ(t),Ω(t))，xmis,i的條件分布為N(ci,σi)，其中

于是，可推得[4,10-11]

(6)

以及

式中：Rij=1----xij為觀測數(shù)據；

Rij=0----xij缺失，i=1,2,…,n，j,k=1,2,…,p。

2)M步。

M步通過最小化Q(μ,Ω|μ(t),Ω(t))來更新參數(shù)估計(μ(t+1),Ω(t+1))，

(8)

式(8)中Ω(t+1)的優(yōu)化問題可直接由GLasso算法求解。將(μ(t+1),Ω(t+1))作為下一次迭代的當前參數(shù)，重復E步和M步直至算法收斂。

m=1,2,…,k，

(9)

其中

3 模擬研究

在完全隨機缺失的機制下，文中對自適應MissGLasso方法和MissGLasso方法在高斯圖模型上的參數(shù)估計和模型選擇進行了模擬比較?，F(xiàn)考慮以下兩個模型：

1)AR(1)，(Ω-1)jk=0.7|j-k|。

2)Ω=B+δI，其中B的對角線元素為0，非對角線元素之間相互獨立，取值為 0.5或0，且P(Bjk=0.5)=0.1。選擇恰當?shù)摩?，使得Ω的條件數(shù)為p。

對于以上兩個模型，設定隨機變量個數(shù)p=100,200,300，對兩個模型與隨機變量個數(shù)p的6種組合，均產生樣本量n=200的50個獨立的數(shù)據集，且在每個數(shù)據集上完全隨機刪除10%、20%、30%的數(shù)據。因此，每個完全觀測的數(shù)據集對應3個具有不同缺失比例的缺失數(shù)據集。為評價算法性能，文中比較以下指標。

Kullback-Leibler損失(kl)、真陽率(tpr)、陽性預測率(ppv)以及馬修斯相關系數(shù)(mcc)。

式中:Ω----真實的協(xié)方差逆陣;

tp----真陽類的個數(shù);

tn----真陰類的個數(shù);

fp----假陽類的個數(shù);

fn----假陰類的個數(shù)。

基于模型1和模型2的不同參數(shù)設定各評價指標的均值和標準差分別見表1和表2。

表1 不同參數(shù)設定下模型1各評價指標的均值和(標準差)

表2 不同參數(shù)設定下模型2各評價指標的均值和(標準差)

由表中可以看出，對于模型1和模型2的所有設定，除了tpr以外，自適應MissGLasso方法均優(yōu)于MissGLasso，并且模型1的tpr僅在缺失比例為30%時略低。對于模型2，自適應MissGLasso方法的mcc略低，但是兩者相差不大。kl為評價協(xié)方差逆陣估計效果的指標，不論模型1還是模型2，自適應MissGLasso的估計性能都更佳。綜上所述，自適應的方法更優(yōu)。

4 實例分析

基于自適應MissGLasso方法，文中對枯草芽孢桿菌生產核黃素(維生素B2)的數(shù)據集進行了圖模型分析。該數(shù)據集由DSM(瑞士)提供，可在文獻[12]的補充材料中下載。數(shù)據為完全觀測，包含71個樣本，4 088個協(xié)變量以及1個響應變量。協(xié)變量為4 088個基因表達水平的對數(shù)；響應變量為核黃素生產率的對數(shù)。為簡單展現(xiàn)自適應MissGLasso方法的圖模型選擇效果，文中僅對經驗方差最大的100個基因，以及測量核黃素生產率的響應變量作圖模型推斷，即數(shù)據集riboflavinv100.csv[12]。該數(shù)據集樣本量n=71，隨機變量個數(shù)p=101。

不同缺失程度下重疊邊數(shù)的箱線圖如圖1所示。

圖1 不同缺失比例下的重疊邊數(shù)

從圖1中可以看出，數(shù)據缺失的比例越高，重疊的邊數(shù)越少。即使數(shù)據集中缺失30%的數(shù)據，自適應MissGLasso方法仍能識別出約30條重疊的邊。

不同缺失比例下的平均稀疏矩陣如圖2所示。

由圖2可以看出，該矩陣的(j,k)元素為50次模擬中協(xié)方差逆陣估計值(j,k)元素的非零頻率。

圖2說明缺失比例越高,選出的邊越少，圖模型越稀疏。

圖2 不同缺失比例下的平均稀疏矩陣

5 結 語

基于懲罰似然框架提出自適應MissGLasso方法以處理缺失數(shù)據情形下協(xié)方差逆陣的估計問題，采用EM算法和GLasso算法進行優(yōu)化求解。該方法可用于數(shù)據完全隨機缺失的情況。通過不同模型下的模擬實驗結果顯示，自適應MissGLasso方法的圖模型選擇性能及協(xié)方差估計結果較MissGLasso方法更優(yōu)。此外，對核黃素生產數(shù)據集的圖模型結構學習實驗同樣驗證了自適應MissGLasso方法具有良好的模型選擇性能。