帶狀結(jié)構(gòu)的高斯圖模型嵌套懲罰估計(jì)

李凡群，張新生

（1.安徽財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030；2.復(fù)旦大學(xué) 管理學(xué)院，上海 200433）

0 引言

無(wú)向圖可以用G=(V，E)來(lái)描述，其中V包含p維隨機(jī)變量X=(X1，…，Xp)的p個(gè)頂點(diǎn)（node），邊集(edge)E={eij|1≤i＜j≤p} 描述了X=(X1，…，Xp)這p個(gè)頂點(diǎn)之間的兩兩條件相依性。當(dāng)且僅當(dāng)Xi和Xj在給定其他變量時(shí)是條件獨(dú)立的，則稱(chēng)Xi和Xj之間不存在邊，即E中不存在元素eij。如果X=(X1，…，Xp)服從高斯分布N(μ，Σ)，則記 Ω=Σ-1為協(xié)方差矩陣的逆矩陣，稱(chēng)為精確矩陣。在高斯分布場(chǎng)合的圖模型中，Xi和Xj之間不存在邊等價(jià)于給定其余變量的條件下，Xi和Xj是條件獨(dú)立的，這又等價(jià)于Ω中元素ωij=0。于是高斯圖模型的恢復(fù)和參數(shù)估計(jì)問(wèn)題就等價(jià)于其協(xié)方差矩陣或其精確矩陣的估計(jì)以及精確矩陣中的零元素的識(shí)別問(wèn)題。

Lauritzer[1]介紹了高斯圖模型的性質(zhì)，并把圖模型選擇和估計(jì)的方法應(yīng)用于高斯圖模型。Meinshausen和Bulmann[2]提出相鄰點(diǎn)選擇的方法。Levina等和Fredman等[3,4]把懲罰似然的方法用于精確矩陣的估計(jì)及圖模型的恢復(fù)。Fan和Peng[5]用自適應(yīng)Lasso和SCAD懲罰選擇圖模型和參數(shù)估計(jì)。Zhang等(2014)[6]提出D-trace損失，在約束條件下，通過(guò)最小化D-trace損失，獲得了稀疏的精確矩陣的估計(jì)。

尤其是多維隨機(jī)變量有序的，比如縱向數(shù)據(jù)、時(shí)間序列、空間數(shù)據(jù)及光譜數(shù)據(jù)等。該類(lèi)數(shù)據(jù)的特點(diǎn)是頂點(diǎn)間的條件相依性與頂點(diǎn)之間的位置有關(guān)，位置越遠(yuǎn)，條件相依性越弱，稱(chēng)該類(lèi)模型為帶狀結(jié)構(gòu)高斯圖模型。對(duì)這類(lèi)數(shù)據(jù)，Smith和 Kohn[7]，Huang等[8]提出對(duì)正定矩陣做修正的Cholesky分解：Ω=(1-Φ)TD-2(1-Φ)，其中Φ是主對(duì)角線(xiàn)上的元素為1的下三角矩陣。D為對(duì)角形矩陣。Bickel等(2008)[9]通過(guò)對(duì)協(xié)方差矩陣的元素設(shè)置門(mén)限以確定帶寬，Levina等(2008)[3]對(duì)修正的Cholesky分解因子Φ提出嵌套的Lasso估計(jì)。Cai等(2011)[10]提出l1懲罰CLME方法，通過(guò)解線(xiàn)性規(guī)劃，獲得協(xié)方差矩陣的估計(jì)。本文利用精確矩陣的諸列(行)元素的系列線(xiàn)性回歸模型的解釋?zhuān)瑢?duì)系列回歸模型的回歸系數(shù)施加嵌套的懲罰，得到精確矩陣的估計(jì)。與有關(guān)文獻(xiàn)相比，該方法更能直接獲得估計(jì)結(jié)果。隨機(jī)模擬結(jié)果表明，對(duì)精確矩陣施加嵌套的懲罰的估計(jì)提高了估計(jì)的精度。

1 精確矩陣的諸列元素的系列線(xiàn)性回歸模型的解釋

令X=(X1，…，Xp)T為一個(gè)p維隨機(jī)變量，且X～N(μ，Σ)。為了表達(dá)方便，本文用Σi，j表示移去Σ的第i行和第j列的的子矩陣，Σi，j表示Σ的第i行，其中第j個(gè)元素被移去，Xi表示移去變量X的第i個(gè)變量后的p-1維的隨機(jī)向量。如果對(duì)X，μ，Σ做如下分塊：

則由多元統(tǒng)計(jì)知，在給定X(2)=x(2)的條件分布為[11]：

特別地，取X(1)=Xi，X(2)=Xi

那么在給定其余變量Xi=xi的條件下，Xi的條件分布為：

如果把Xi看成響應(yīng)變量，其余變量Xi看成協(xié)變量，考慮如下的線(xiàn)性回歸模型：

此時(shí)令β(j)=(β1，…，βj-1，βj+1…，βn)T，?i與Xi相互獨(dú)立，并且在均方意義下：

另一方面，若令Ω為多元正態(tài)分布X～N(μ，Σ)的精確矩陣，則由ΣΩ=I，有：

對(duì)每一個(gè)i(i=1，2，...，p)，得到精確矩陣Ω的諸對(duì)角線(xiàn)上的元素ωii以及其第i列上其余元素Ωi，i與回歸問(wèn)題（1）中回歸系數(shù)β(j)及殘差方差之間有如下關(guān)系：

這表明精確矩陣Ω的第i列的元素Ωi，i與回歸系數(shù)成比例，于是Ωi，i中的零元素與β(j)中的零元素相對(duì)應(yīng)，反之亦然。從而高斯圖模型的參數(shù)估計(jì)與模型選擇可以轉(zhuǎn)化為對(duì)回歸問(wèn)題（1）的模型選擇及對(duì)殘差的方差的估計(jì)。基于估計(jì)

2 精確矩陣嵌套的l1懲罰估計(jì)

2.1 嵌套懲罰函數(shù)

根據(jù)前文的分析知，精確矩陣Ω中每行元素有相應(yīng)的回歸問(wèn)題中回歸系數(shù)的解釋。并且模型恢復(fù)與參數(shù)估計(jì)的問(wèn)題轉(zhuǎn)化為對(duì)回歸系數(shù)β(j)和殘差方差的選擇和估計(jì)問(wèn)題。

基于模型（1）的β(j)和的負(fù)的對(duì)數(shù)似然函數(shù)為：

(j)懲罰似然的估計(jì)為：

Levina等(2008)[3]基于改進(jìn)的Cholesky分解，對(duì)cholsesky分解因子每行元素提出嵌套的懲罰，得到嵌套的分解因子估計(jì)。

下面本文基于帶狀模型對(duì)應(yīng)的回歸模型，對(duì)其系數(shù)β(j))提出如下的嵌套懲罰函數(shù)：

該懲罰函數(shù)相當(dāng)于加權(quán)的l1懲罰函數(shù)。在該懲罰下，一方面，當(dāng) k＜ j 時(shí)，若0。這意味著精確矩陣Ω的第j行元素有估計(jì)當(dāng)l＞j時(shí)，若這意味著精確矩陣Ω的第j行元素有估計(jì)于是通過(guò)該懲罰，得到的圖模型估計(jì)保證了其帶狀結(jié)構(gòu)特性。

考慮到 |βj-1|，|βj+1|的分母為1，故對(duì) |βj-1|，|βj+1|兩個(gè)參數(shù)的懲罰在尺度上和其余參數(shù)不一致，所以對(duì)P1(|β(j)|)做如下兩種改進(jìn)：

2.2 局部平方近似以及牛頓迭代算法

本文利用迭代算法，通過(guò)最小化：

得到β(j)和估計(jì)。

第一步：給定β(j)的當(dāng)前估計(jì)值,計(jì)算：

對(duì)每一個(gè)j=1，2，…，p重復(fù)以上兩步，直到滿(mǎn)足一定的準(zhǔn)則，從而得到諸β(j)和的估計(jì)，根據(jù)前文的分析，得到精確矩陣的估計(jì)。由于函數(shù)不是關(guān)于β(j)的凸函數(shù)，為了克服計(jì)算上的困難，本文對(duì)諸進(jìn)行局部平方近似，得到了的顯式表達(dá)式。

一方面，對(duì)于懲罰函數(shù)的分子部分，令：

則上述三個(gè)懲罰函數(shù)經(jīng)過(guò)局部平方近后得到如下形式：

其中C1為常數(shù)為對(duì)角形矩陣:

其中C2為常數(shù)，為對(duì)角形矩陣:

其中C3為常數(shù)，為對(duì)角形矩陣:

3 隨機(jī)模擬

下面通過(guò)幾個(gè)例子，分別計(jì)算估計(jì)的損失，比較不同方法的估計(jì)的精度。本文用K-L損失來(lái)度量估計(jì)的精度：把本文的方法得到的估計(jì)Modified nested Lasso簡(jiǎn)記為MNL，并把結(jié)果同以下方法的估計(jì)進(jìn)行比較，包括基于Cholesky分解的nested-Lasso(CNL)[3]，g-Lasso[4]，以及由經(jīng)驗(yàn)協(xié)方差矩陣得到的精確矩陣的估計(jì)記為Sample?？紤]p=30，50，100，200的情況，對(duì)于每一個(gè)例子中的每一個(gè)p，從真實(shí)的模型中產(chǎn)生50組樣本數(shù)據(jù)，每組樣本數(shù)據(jù)各包含200個(gè)觀(guān)測(cè)值。分別利用BIC和CV準(zhǔn)則，選取調(diào)節(jié)參數(shù)λ。

在進(jìn)行迭代過(guò)程中，當(dāng)參數(shù)的估計(jì)的絕對(duì)值|βk|＜10-5時(shí)，令 |βk|=0。由于模型1是稀疏的，所以本文報(bào)告了BIC準(zhǔn)則和CV準(zhǔn)則下的估計(jì)結(jié)果。對(duì)于每一個(gè)模型，由式（2）至式（4），本文分別計(jì)算了三種迭代估計(jì)的結(jié)果，列出其中估計(jì)的最小損失以及損失標(biāo)準(zhǔn)差。所有結(jié)果列于表1和表2所示。表1和表2表明，基于精確矩陣元素的回歸模型的改進(jìn)的嵌套懲罰估計(jì)一致優(yōu)于所列的其他方法。

表1 K-L損失下不同的估計(jì)方法對(duì)模型1的估計(jì)結(jié)果

表2 K-L 損失下不同的估計(jì)方法對(duì)模型2 的估計(jì)結(jié)果

（1）模型1(AR(3))：精確矩陣Ω=(ωij),其中ωij=0.5|i-j|，|i-j|≤3，i≠j,以及ωii=1,其他場(chǎng)合ωij=0。

（2）模型2（完全模型）：精確矩陣 Ω=(ωij),其中ωij=0.5|i-j|。

4 結(jié)束語(yǔ)

本文在具有帶狀結(jié)構(gòu)的高斯圖模型場(chǎng)合下，利用高斯圖模型的精確矩陣的諸元素可以利用系列線(xiàn)性回歸模型進(jìn)行解釋的性質(zhì)，對(duì)諸線(xiàn)性回歸系數(shù)施加嵌套的l1懲罰，得到精確矩陣的改進(jìn)的嵌套Lasso估計(jì)。通過(guò)局部平方近似及牛頓迭代算法，能很方便得到估計(jì)結(jié)果。與基于改進(jìn)的Cholesky分解的嵌套Lasso估計(jì)相比，本文的方法能更直接得到精確矩陣的估計(jì)，減少了估計(jì)損失。數(shù)值模擬結(jié)果表明精確矩陣的嵌套Lasso估計(jì)提高了估計(jì)速度和估計(jì)精度。