貝葉斯圖彈網模型及其推廣

葉贇鑫，傅德印，2

（1.蘭州財經大學 統計學院，蘭州 730020；2.中國勞動關系學院，北京 100048）

0 引言

隨著語音處理、圖像處理等領域的不斷發展，高維復雜數據越來越常見，變量之間的關系也逐漸變得復雜，圖模型提供了解決這些問題的一般方法[1]，其主要通過圖表的形式來增強對復雜概率模型的理解與學習；同時提供了一種簡單的方法來可視化概率模型的結構，通過檢查圖可以輕松了解到模型的各項屬性[2]。而高斯圖模型是一種特殊的圖模型[3,4]。對精度矩陣的估計和模型恢復的研究通過在圖模型中施加LASSO懲罰來達到目的[5],這種方法也被稱為圖LASSO，通過最大化懲罰對數似然實現參數估計。李凡群和張新生（2018）[6]對各線性回歸系數施加嵌套的L1懲罰,得到精度矩陣改進的嵌套LASSO估計；李凡群和張孔生（2021）[7]利用自適應LASSO 估計的貝葉斯解釋,構造了懲罰函數的分層貝葉斯先驗模型；鄭倩貞等（2021）[8]考慮了在潛變量高斯圖模型下的模型選擇問題，介紹了高斯圖模型及在潛變量高斯圖模型下的LVglasso方法。在此基礎上，Wang（2012）[9]提出了貝葉斯圖LASSO，并將其推廣到貝葉斯自適應圖LASSO，與圖LASSO[5]、自適應圖LASSO、SCAD[10]等相比，整體上表現出明顯的優勢。然而LASSO 方法存在著一定的缺陷。Zou 和Hastie（2005）[11]在LASSO 的基礎上增加了L2范數，提出了彈網（Elastic Net）模型以彌補LASSO的不足。其中，L1范數提供了稀疏可解釋性，而L2范數提高了方法的靈敏度[12]。因此，在高斯圖模型中加入彈網懲罰也是一個可行的方案。在研究圖模型時，為了保持類內關系并減少類間關系，Liu等（2016）[13]提出了約束彈網正則化模型；Bernardini等（2021）[14]、Kovács 等（2021）[15]在高斯圖模型的背景下將圖LASSO模型推廣至圖彈網模型。

考慮到彈網模型的“雙縮問題”以及貝葉斯方法在統計推斷上的簡潔性，Li 和Lin（2010）[16]提出了貝葉斯彈網模型。為了進一步提升高斯圖模型對精度矩陣的估計以及模型恢復的能力，本文考慮在圖彈網的基礎上將貝葉斯彈網模型引入高斯圖模型。在貝葉斯圖LASSO和貝葉斯自適應圖LASSO[9]的基礎上提出貝葉斯圖彈網（Bayesian Graphical Elastic Net，BGEN）模型并推廣至貝葉斯自適應圖彈網（Bayesian Adaptive Graphical Elastic Net，BAGEN）模型；同時，介紹了圖彈網模型的完全貝葉斯處理以及彈網懲罰下的分塊Gibbs（Block Gibbs）采樣方法[9]。

1 方法介紹

1.1 貝葉斯圖LASSO模型

在高斯圖模型的研究中，Yuan 和Lin（2007）[17]通過LASSO 懲罰提升了其模型選擇與參數估計的能力。這種方法也被稱為圖LASSO，通過最大化懲罰對數似然實現參數估計：

其中，Θ=Σ-1為精度矩陣，S為經驗協方差矩陣。隨后，Wang（2012）[9]在貝葉斯框架下進行了推廣，通過其提出的分塊Gibbs方法實現了圖LASSO的完全貝葉斯處理，提出了貝葉斯圖LASSO模型和貝葉斯自適應圖LASSO模型。相較于頻率派和其他貝葉斯方法，貝葉斯圖LASSO有更好的矩陣估計和圖形結構學習性能。Wang（2012）[9]對精度矩陣Θ的對角元素和非對角元素分別賦予了獨立的指數先驗和雙指數先驗，因此貝葉斯圖LASSO 模型相當于以下模型的最大后驗估計：

其中，DE(x|λ)是形式為p(x)=λ/2 exp(-λ|x|)的雙指數密度函數，Exp(x|λ)是形式為p(x)=λexp(-λx)Ix＞0的指數密度函數，N(y|，Σ)為多元正態分布，C是與懲罰參數λ無關的歸一化常數，I為指示函數。同時，為了使式（2）能夠通過Gibbs方法進行采樣，需將其轉換成正態分布的尺度混合分布。因此，貝葉斯圖LASSO模型的分層模型為：

貝葉斯框架下傳統的MCMC 方法難以對精度矩陣Θ進行有效估計，而將精度矩陣Θ的條件分布適當地參數化之后，分塊Gibbs 采樣將是一個十分有效的采樣器[9]。根據式（3），貝葉斯圖LASSO在使用分塊Gibbs采樣時，目標分布如下：

其中，精度均值Θ被約束為正定矩陣。同時，當使用分塊Gibbs采樣時，需要對精度矩陣Θ進行分塊處理。令T=(τij)為p×p的對稱矩陣，對角線元素為0。得到分塊矩陣Θ、S和T：

其中，θi=(θi，1，θi，2，…，θi，i-1，θi，i+1，…，θi，p)T為精度矩陣Θ第i列的所有非對角元素所組成的向量，θii為精度矩陣Θ對角線上的第i個元素，Θii為除去精度矩陣Θ的第i行和第i列元素后得到的(p-1) ×(p-1) 矩陣。同理可得，sii和0分別為矩陣S和矩陣T對角線上的第i個元素，si和τi分別為矩陣S和矩陣T第i列的所有非對角元素所組成的向量，Sii和Tii分別為去除矩陣S和矩陣T的第i行和第i列元素后得到的(p-1) ×(p-1) 矩陣。

根據式（5）的3個分塊矩陣，Θ第i列的條件分布為：

其中，Dτ=diag(τi) 。在式（6）的基礎上，Wang（2012）[9]將變量進行轉換，令β=θi,，則有：

根據式（7），(γ，β)的條件分布為：

對于懲罰參數λ，Wang（2012）[9]給定了一個伽馬先驗分布λ～Gamma(r，s)，得到其條件后驗分布為：

其中，r和s為給定的常數。最后，Wang（2012）[9]通過對精度矩陣的所有元素賦予不同的懲罰參數，將貝葉斯圖LASSO模型拓展至貝葉斯自適應圖LASSO模型，進一步提升了模型的穩健性。貝葉斯自適應LASSO的后驗形式為：

1.2 貝葉斯圖彈網模型

1.2.1 圖彈網先驗

對于圖彈網模型，本文考慮對精度矩陣Θ的所有非對角元素施加一個形式為p(x)=λ1/2 exp(-λ1|x|-λ2x2)的彈網先驗[16]。因此，本文所考慮的圖彈網具有以下形式的最大后驗估計：

其中，EN(x|λ1，λ2)是形式為p(x)=λ1/2 exp(-λ1|x|-λ2x2)的彈網先驗。為了將圖彈網模型推廣至貝葉斯框架并使其能夠通過Gibbs采樣方法對模型進行采樣，需要將式（13）中的彈網先驗進行轉換。通常需要將EN(x|λ1，λ2)替換成正態分布的尺度混合[16]。然而對θij使用正態分布的尺度混合分布時，由于精度矩陣必須是正定矩陣，因此θij將不再獨立于給定的尺度參數[9]。精度矩陣Θ的上三角元素所組成的向量與潛在尺度參數具有以下形式的條件概率密度函數：

其中，規范化常數Cτ依賴于τ且沒有解析表達式。同時，選取p(Θ|τ)和p(τ|λ1，λ2)作為先驗，由于：

雖然規范化常數Cτ沒有解析表達式，但根據式（14）和式（15）進行分塊Gibbs采樣時可以將其抵消。因此，θij的邊際分布依舊遵循式（13）的形式。

1.2.2 分塊Gibbs采樣

在貝葉斯框架下，傳統的MCMC方法很難對精度矩陣Θ進行有效估計[18]。而將精度矩陣Θ的條件分布適當地參數化之后，分塊Gibbs 采樣將是一個十分有效的采樣器[9]。因此，貝葉斯圖彈網在使用分塊Gibbs 進行采樣時具有以下目標分布：

其中，矩陣S=YTY/n，Y=(y1，y2，…，yn)是由p個變量和n個樣本量組成的n×p矩陣。同時，精度矩陣Θ為正定矩陣，潛在尺度參數τ沒有限制。為了能夠使用分塊Gibbs 采樣對精度矩陣進行估計，需要對式（17）中的矩陣進行處理。令T=(τij為p×p的對稱矩陣，其對角線上的元素全部等于0。矩陣Θ、S和T的分塊處理如下：

根據式（17）的分布形式，Θ第i列的條件分布具有以下形式：

其中，Hτ=2λ2diag[τi/(τi-1)]。為了方便計算，需要對變量進行轉變。令η=θi，。則：

根據式（20）可以得到η和κ的條件分布為：

對于潛在尺度參數τ，其迭代更新與貝葉斯彈網[16]類似。1/ (τij-1) 的條件后驗分布獨立且服從逆高斯分布：

綜上所述，根據式（17）至式（23），進一步將貝葉斯圖彈網模型的分塊Gibbs 采樣總結為算法1（分塊Gibbs 算法）。過程如下：

輸入：對于第t次迭代，給定當前值Θ(t-1)∈M+，τ(t-1)。

步驟1：fori=1，2，…，p，進行下面的迭代，直到收斂：

根據式（18），對矩陣Θ、S和T進行分塊處理；

根據式（21）和式（22），采樣得到η和κ；

步驟2：當i＜j時，從采樣得到aij，更新=(1/aij)+1。

輸出：Θ(t)∈M+，τ(t)。

1.2.3 懲罰參數的選擇

對于圖彈網中懲罰參數λ1和λ2的選擇，常用的方法是最大化數據邊際似然的經驗貝葉斯估計[16]，或者使用交叉驗證對懲罰參數λ1和λ2進行選擇[11,14,15]。但是對于貝葉斯框架下的圖彈網模型來說，使用以上兩種方法并不是完全貝葉斯方法。而在貝葉斯圖LASSO模型中已經實現了通過完全的貝葉斯分析方法對懲罰參數進行估計[9]。根據式（13），規范化常數C具有以下積分形式:

可以注意到，式（24）中規范化常數C并沒有解析表達式且依賴于懲罰參數λ1和λ2。為了能夠從滿條件后驗分布中對λ1和λ2進行采樣,可以考慮用與式（14）和式（15）類似的方式將規范化常數C抵消[9]。因此，對懲罰參數λ1和λ2給定以下先驗：

對于式（25）中的超參數，設置αλ1=1，αλ2=0.01*p和βλ1=βλ2=0.01，其中，p為精度矩陣的維度。根據式（25）及式（13）可以得到條件后驗為：

1.2.4 模型檢驗方法

圖彈網模型的主要目標是對數據的協方差矩陣以及精度矩陣進行參數估計和變量選擇。而其估計的結果可以通過Stein 損失函數進行評估。給定協方差矩陣Σ=Θ-1，Stein損失函數可以表示為：

根據蒙特卡洛采樣得到的Θ即可估計。另外，對于圖形結構學習性能的評估，可以使用馬修斯相關系數（Matthews Correlation Coefficients，MCC）[9,10]：

其中，TP、TN、FP和FN分別為真陽性、真陰性、假陽性和假陰性的數量。由于同時考慮了TP、TN、FP和FN，因此MCC 通常被認為是一種平衡的方法。對于MCC，數值越大，說明分類的效果越好，即矩陣的還原效果越好。同時根據Fan等（2009）[10]的研究，對于需要估計的參數θij，若＜10-3，則令{θij=0}。

2 模型拓展

在常見的變量選擇方法中，自適應LASSO 模型具有SCAD 的Oracle 性質[20]，但是對于高維數據的處理缺乏穩定性；彈網模型能夠很好地處理高維情況下的共線性問題[11]，但不滿足Oracle性質。而自適應彈網模型結合了上述模型的優點，在處理共線問題方面優于其他算法，大大提高了有限樣本的性能[21]。本文考慮給予精度矩陣中的每個非對角元素θij不同的懲罰參數λ1和λ2。將貝葉斯圖彈網模型推廣至貝葉斯自適應彈網模型。

根據式（13）和式（25），貝葉斯圖彈網模型具有和式（11）類似的處理規范化常數的方式。同樣，自適應圖彈網模型給定精度矩陣Θ的非對角元素一個形式為p(x)=λ1/2 exp(-λ1|x|-λ2x2)的彈網先驗EN(x|λ1，λ2)：

其中，規范化常數C{λ1，λ2，λ}在計算條件后驗時被抵消了。同時，在式（31）中，精度矩陣Θ的每個非對角元素都被給予了不同的懲罰參數。因此，對于任意的i＜j，懲罰參數λ1，ij和λ2，ij服從不同形式的伽馬分布：

根據伽馬分布的性質，θij較小時將會得到一個較大的懲罰參數λ1和λ2，反之亦然。但是，為了降低先驗分布的超參數對θij的影響，提升懲罰參數λ1和λ2對θij的適應性，式（32）中超參數必須足夠小[9]。因此將自適應圖彈網的超參數設定為α1=0.01，α2=0.005*p，β1=β2=10-4，其中，p為精度矩陣的維度。

3 數值模擬

為了驗證本文方法的有效性，根據Wang（2012）[9]的研究，本文考慮了6個不同的模型。對于所考慮的6個不同的模型，分別設置了不同的協方差矩陣和精度矩陣，并通過這6 個模型檢驗所提模型在不同情形下的表現。本文所使用的所有模型相關參數均服從正態分布：

下面給出本文考慮的6 個模型，這些模型參照了Wang（2012）[9]的假設，只設定了協方差矩陣或者精度矩陣：

模型1：令協方差矩陣Σ中所有元素σij=0.7|i-j|。

模型2：令精度矩陣Θ的對角線元素θii=1，非對角線元素θi，i-1=θi-1，i=0.5 和θi，i-2=θi-2，i=0.25。

模型3：令協方差矩陣Σ的對角線元素σii=1。當1 ≤i≠j≤p/2 或者p/2 ≤i≠j≤p時，非對角線元素σij=0.5。

模型4：令精度矩陣Θ的對角線元素θii=1，第一行與第一列元素θi，1=θ1，i=0.1，其余元素全都為0。

模型5：令精度矩陣Θ的對角線元素θii=2，非對角線元素θi，i-1=θi-1，i=0.5 與θp，1=θ1，p=0.9。

模型6：令精度矩陣Θ的對角線元素θii=2，所有非對角線元素θij=1。

對于以上6個模型，每個模型均考慮生成一組維度為p=30、樣本量為n=100 的數據，一組維度為p=100、樣本量為n=500 的數據以及一組維度為p=100、樣本量為n=80 的數據。對每組數據進行50 次模擬試驗并計算每次模擬的Stein損失以及對應的MCC值。本文一共檢驗對比了4 種不同的帶正則化的圖模型。給出50 次模擬試驗的Stein損失的中位數和對應標準差結果（如表1所示）。

從表1可知，當p＜n時，兩種自適應模型整體上要比非自適應模型表現更優。從矩陣維度的角度來看，在維度較低時，圖彈網模型和圖LASSO模型結果相似，自適應圖彈網模型和自適應圖LASSO 模型結果相似，但是兩種圖彈網模型的Stein損失結果的標準差相對兩種圖LASSO模型要小一些；在維度較高的情況下，貝葉斯自適應圖彈網模型的性能是最好的，貝葉斯圖彈網模型的性能優于貝葉斯圖LASSO模型。從模型的角度來看，針對模型3，4種模型其實都有著相似的性能；針對模型6，兩種非自適應模型則要優于兩種自適應模型；針對其他4 個模型，兩種自適應模型則是要明顯優于另外兩種非自適應模型。當p＞n時，4種模型在對模型5進行分析時均有極大可能導致協方差矩陣非正定，從而使得模型失效，因此本文并未得到模型5的結果。通過對比其余5個模型下的模擬結果可以發現，貝葉斯自適應圖彈網模型整體上要優于其余模型。因此，對于圖模型的參數估計與變量選擇，貝葉斯自適應圖彈網模型幾乎是整體最優的，尤其是在矩陣維度較高的情況下。帶彈網懲罰的圖模型的參數估計與變量選擇結果的穩定性更強；同時，自適應的模型也要比非自適應的模型更優。針對模型6這種全模型來說，非自適應的模型則要更好。模型4 在p=100 的情況下，貝葉斯圖LASSO模型的結果非常差。

對于MCC，得到的數值越大說明模型的還原效果越好。表2給出了50次模擬的MCC的均值與標準差。

表2 MCC結果

由于模型6 為全模型，MCC 結果全都等于零，因此，表2 只給出了模型1 至模型5 的MCC 結果。根據表2的結果可以看出，對于所有模型，MCC結果的方差都非常小，且從整體上來說，貝葉斯自適應圖彈網模型的表現依舊是最好的。同樣，對于圖模型的分類效果，貝葉斯自適應圖彈網和貝葉斯自適應圖LASSO模型也要優于另外兩種非自適應模型。另外，模型3 是所有模型中分類效果最好的模型，且4 種對比模型并沒有太大的區別。

4 實例分析

本文通過幾組實例數據，進一步驗證貝葉斯圖彈網模型在參數估計與模型選擇方面的優越性和對真實數據的適應性。使用以下3組數據來對模型進行檢驗：第1組數據是由Sachs 等（2005）[22]收集的流式細胞術數據，其中包含了n=7466 個細胞和p=11種蛋白。第2組數據是通過收集2020 年1 月至2021 年6 月上證指數與深證綜指的p=35 只消費股與銀行股，共得到n=384 個交易日的日回報數據。第3 組數據為n=60 個來自美國猶他州的具有北歐與西歐血統的無關個體的基因數據[23]，包含了p=100個變量。

為了比較模型性能，考慮對每組數據進行5折交叉驗證，并通過Stein 損失與MCC 進行度量。Stein 損失函數具有以下形式：

表3 實例分析Stein損失結果

表4 實例分析MCC結果（單位：%）

由表3 與表4 可知：首先，針對第1 組數據，當數據的樣本量遠大于變量時(n?p)，本文所比較的四種模型的結果非常接近；且Stein損失和MCC的結果表明，這4種模型均有十分良好的性能。因此可以認為在這種情況下4種模型都適用。其次，針對第2 組數據，當數據的變量與樣本量比較接近時(n＞p)，貝葉斯自適應圖彈網和貝葉斯自適應圖LASSO 模型是要優于另外兩種模型。而4 種模型中性能最優的是貝葉斯自適應圖彈網模型。最后，針對第3 組數據，當數據的變量大于樣本量時(n＜p)，兩種帶彈網懲罰的模型要明顯優于另外兩種帶LASSO懲罰的模型。

綜上所述，在上述4 種模型中，貝葉斯自適應圖彈網模型的適應性最強，對于3 組數據均表現出了良好的性能。此外，貝葉斯圖彈網和貝葉斯自適應圖LASSO 模型在特定類型的數據下也表現出了較好的性能。

5 結束語

本文在貝葉斯圖LASSO 模型的基礎上，使用彈網懲罰來代替LASSO懲罰，提出了貝葉斯圖彈網模型，并將貝葉斯圖彈網拓展至貝葉斯自適應圖彈網；同時，通過分塊Gibbs采樣對精度矩陣進行估計。通過數值模擬和實例分析的結果可知，貝葉斯圖彈網模型在維度較高的情況下對模型的協方差矩陣與精度矩陣的識別能力要優于貝葉斯圖LASSO模型；同時，貝葉斯自適應圖彈網模型在參數估計和模型選擇方面有著良好的性能，在對實例數據進行分析時也有著很強的適應性。