基于自助法的高斯貝葉斯網(wǎng)結(jié)構(gòu)學(xué)習(xí)

徐平峰, 王樹達(dá), 尚來旭, 楊 哲

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

0 引 言

貝葉斯網(wǎng)是概率論和圖論的有機(jī)結(jié)合，其圖結(jié)構(gòu)是一個有向無圈圖(DAG)，可清晰直觀地表示變量間的條件獨立關(guān)系、因果關(guān)系和變量間的依賴關(guān)系。在貝葉斯網(wǎng)領(lǐng)域中，結(jié)構(gòu)學(xué)習(xí)是一類非常重要的研究問題，即從觀測數(shù)據(jù)出發(fā)探求變量間的條件獨立關(guān)系。

結(jié)構(gòu)學(xué)習(xí)方法大致分為三類:

1)基于約束的方法。這類方法依賴于條件獨立關(guān)系檢驗，其典型代表是PC算法[1]。該算法首先進(jìn)行低階條件獨立性檢驗，然后逐步進(jìn)行高階條件獨立性檢驗，排除不可能的邊，建立有向無圈圖的骨架，確定v結(jié)構(gòu)，最后盡可能地為邊添加方向，找到一個等價類。

2)基于評分搜索的方法。該方法為每個DAG賦一個評分，然后在搜索空間內(nèi)搜索一個得分最低的DAG或等價類。例如,Chickering[2]提出了貪婪等價類搜索(GES)方法，通過逐步向前和向后搜索得到得分最小的等價類；在誤差等方差的條件下，Peters 等[3]利用貪婪DAG搜索方法(GDS)求BIC評分最小的高斯貝葉斯網(wǎng)。若先利用l1懲罰似然學(xué)習(xí)高斯圖模型的無向圖，然后在貪婪搜索時保證DAG的骨架是該無向圖的邊，這時搜索DAG的范圍將變小，該方法稱為GDSM。

3)混合方法。分兩個階段：

a)限制階段，先進(jìn)行一些條件獨立性檢驗，利用檢驗結(jié)果排除大量的邊，以減小候選DAG的空間；

b)得分搜索階段，在這些條件獨立性關(guān)系的約束下，搜索評分最低的DAG。例如，Tsamardinos等[4]提出MMHC算法，先采用基于約束的方法找到一個頂點可能的父頂點和子頂點，然后搜索BDe得分最小的DAG。

近些年，基于自助法(Bootstrap)的模型聚合方法已經(jīng)應(yīng)用到很多領(lǐng)域，如高維數(shù)據(jù)的變量選擇(Bach[5]，Wang等[6])和高斯圖模型的結(jié)構(gòu)學(xué)習(xí)(Meinshausen等[1]，Li[8]);Elidan[9]利用Bootstrap數(shù)據(jù)集計算穩(wěn)健的模型選擇準(zhǔn)則。對于每個數(shù)據(jù)集和一個給定的DAG，求此DAG參數(shù)的極大似然估計，進(jìn)而利用聚合的信息準(zhǔn)則來衡量G;Wang 等[10]提出了一種基于Bootstrap的DAG結(jié)構(gòu)學(xué)習(xí)方法，稱為DAGBag。該方法基于數(shù)據(jù)的Bootstrap重采樣學(xué)習(xí)DAG集合，以結(jié)構(gòu)海明距離(SHD)為度量，通過最小化到整個DAG集合的總距離來導(dǎo)出聚合DAG。

文中提出一種新的基于自助法的結(jié)構(gòu)學(xué)習(xí)算法，并將其應(yīng)用于高斯貝葉斯網(wǎng)。算法分為三個步驟：

1)對觀測樣本重抽樣得到B組Bootstrap樣本；

2)利用結(jié)構(gòu)學(xué)習(xí)算法,如PC算法、GES算法等，得到B個DAG，并求得B個DAG對應(yīng)貝葉斯網(wǎng)的極大似然估計；

3)以B個DAG中出現(xiàn)的邊作為候選邊搜索一個DAG，使其對應(yīng)的貝葉斯網(wǎng)的極大似然估計與B個極大似然估計的懲罰KL距離的平均值最小。

文中算法與Wang 等[11]最主要的區(qū)別在于，利用懲罰的KL距離作為聚合準(zhǔn)則。KL距離比SHD有更好的性質(zhì)，例如對于兩個等價的DAG，它們對應(yīng)模型的極大似然估計是相同的，因而這兩個極大似然估計到B個極大似然估計的懲罰KL距離相同，但是SHD卻不具備這一性質(zhì)。

1 高斯貝葉斯網(wǎng)

貝葉斯網(wǎng)結(jié)構(gòu)由有向無圈圖G(V,E)表示。頂點集V={1,2，…,p}表示隨機(jī)向量X=(X1,X2，…,Xp)。邊集E={(i,j)∈V×V:i→j}, 其中i→j表示圖G中的有向邊。給定圖G的結(jié)構(gòu), 隨機(jī)向量X=(X1,X2,…,Xp)的聯(lián)合概率密度函數(shù)可以分解為

給定一個聯(lián)合分布，就可能存在多種分解形式，從而得到不同的DAG。所有蘊含相同條件獨立關(guān)系的DAG構(gòu)成一個等價類，同一等價類中的DAG具有相同的骨架和v結(jié)構(gòu)。每個等價類都可以由完全部分有向無圈圖(CPDAG)來表示。

當(dāng)聯(lián)合分布P為多元高斯分布時，高斯貝葉斯網(wǎng)可根據(jù)密度函數(shù)的分解寫成結(jié)構(gòu)方程模型的形式。

j=1,2,…,p，

X=v+ATX+E,

其中

v=(v1,v2,…,vp)T，

E=(ε1,ε2,…,εp)～N(0,D)，

A=(aij)是鄰接矩陣,Aii=0。以I表示p×p維單位矩陣，通過簡單的矩陣變換有

X=(I-AT)-1v+(I-AT)-1E

由E～N(0,D)及正態(tài)分布的性質(zhì)，可得X～N(μ,Σ)，其中

μ=(I-AT)-1ν,

Σ=(I-AT)-1D(I-A)-1,

這里

|I-A|=1，

考慮在給定樣本的情況下，如何估計μ、A和D。

令x(1),x(2),…,x(n)為多元高斯分布N(μ,Σ)的n個觀測值，其中

x(i)=(xi1,xi2,…,xip)T。

樣本均值

樣本協(xié)方差陣

則有對數(shù)似然函數(shù)

這里

Σ=(I-AT)-1D(I-A)-1，

于是有

式中：Q=(I-A)D-1(I-A)T。

2 基于Bootstrap和KL距離的結(jié)構(gòu)學(xué)習(xí)

許多貝葉斯網(wǎng)結(jié)構(gòu)學(xué)習(xí)算法不具有穩(wěn)健性，即使是較小的數(shù)據(jù)擾動，都可能導(dǎo)致得到的DAG有較大變化。為此，Wang 等[10]提出了自助聚合(Bootstrap aggregating)方法，以得到更穩(wěn)定的圖結(jié)構(gòu)。這種方法可以在一定程度上避免過擬合，并減小學(xué)習(xí)算法的假邊發(fā)現(xiàn)率。它們的算法稱為DAGBag，詳細(xì)過程見表1。

表1 算法1. DAGBag

上述DAGBag算法以DAG之間的結(jié)構(gòu)海明距離(SHD)為度量，旨在尋找一個與集合Ge:={G1,G2,…,GB}中所有DAG的平均距離最小的有向無圈圖G*。但是這種采用SHD為度量的方法不足以區(qū)分每一個DAG，原因在于：等價的DAG可能得到不同的scored(G:Ge)，而不同等價類中的DAG也可能得到相同的scored(G:Ge)。

為避免這個問題，以Kullback-Leibler(KL)[11]距離作為度量。兩個概率密度函數(shù)f(x)和g(x)之間的KL距離定義為

對于正態(tài)分布N(μ1,Σ1),N(μ2,Σ2)，其KL距離為

當(dāng)μ1=μ2時，有

從而有

式中：C0----與樣本協(xié)方差陣S有關(guān)的量。

表2 算法2. 基于KL距離的Bootstrap學(xué)習(xí)DAG

選取KL距離作為度量，是因為：

1)可以更精確地度量分布間的距離。

以及正定矩陣Σ，有

其中

式中：n----樣本量;

df----模型中自由參數(shù)的個數(shù)，也就是鄰接矩陣A中的非零元素個數(shù)。

表3 算法3. 基于Bootstrap的懲罰KL距離學(xué)習(xí)DAG

在算法3中，有以下幾點需要說明。首先，針對算法c)中的優(yōu)化問題

采用類似于GDS的搜索策略。由于DAG數(shù)量隨變量數(shù)呈超指數(shù)增長，一般無法遍歷搜索,以求得上式的最小值，因此,借鑒基于道義圖的方法GDSM思想，利用已得到的Ge:={G1,G2,…,GB}信息，以縮小DAG骨架的搜索范圍。具體對于任意一對頂點i和j，計算有向邊i→j在Ge中出現(xiàn)的頻率fij。給定一個閾值f，令骨架邊集候選集合為{i→j|fij≥f}。也就是說，只有頻率大于閾值的邊才可能作為DAG中的邊。閾值f根據(jù)不同的結(jié)構(gòu)學(xué)習(xí)算法L采用不同的值。文中模擬中，如果用GES算法作為L，則選的閾值為f≤0.70。若PC算法作為L，一般f≤0.30。

進(jìn)行貪婪DAG搜索時，考慮兩種初值的方式：

1)根據(jù)fij對所有的有向邊i→j由大到小排序，然后從空的DAG(即沒有邊的DAG)出發(fā)，按邊的頻率由大到小依次向DAG中添加邊，使DAG中的邊盡量與{i→j|fij≥f}一致，并且保證得到的圖仍然是一個有向無圈圖，該初值稱為最多邊DAG；

2)從Ge:={G1,G2,…,GB}中隨機(jī)選一個圖，該初值稱為隨機(jī)DAG。

為了避免貪婪搜索陷入局部極值，將重復(fù)5次貪婪搜索過程，選擇一個得分最好的結(jié)果作為最終結(jié)果。

3 模擬實驗

通過模擬來比較BPKL算法、GDS算法、基于道義圖的GDSM算法、PC算法和GES算法。

在BPKL算法中，a)重抽樣次數(shù)設(shè)置為B=200;b)結(jié)構(gòu)學(xué)習(xí)算法L采用GES算法和PC算法;c)搜索過程中，骨架候選邊頻率的閾值，對于GES算法，選為f=10%,30%,50%,70%,90%，而對于PC算法，選為f=1%,5%,10%,30%。搜索的初值采用上節(jié)所述的兩種方式。如果結(jié)構(gòu)學(xué)習(xí)算法采用GES，f=10%，初值為最多邊DAG，則算法記為gesb10; 如果初值選為隨機(jī)DAG，則算法記為gesbR10;如果結(jié)構(gòu)學(xué)習(xí)算法采用PC，f=5%，初值為最多邊DAG，則算法記為pcb5; 如果初值選為隨機(jī)DAG，則算法記為pcbR5。以此類推。

基于該模型產(chǎn)生樣本量為n的模擬數(shù)據(jù)，并采用以上各種方法進(jìn)行結(jié)構(gòu)學(xué)習(xí)。其中BPKL方法和GDS、GDSM是在DAG的空間中搜索的，因而得到一個DAG；而GES和PC將得到CPDAG，用于表示馬爾可夫等價。為了比較公平性，將前三種方法得到的DAG轉(zhuǎn)換成CPDAG。從而將它們的CPDAG求與真實模型的CPDAG之間的結(jié)構(gòu)海明距離。所謂結(jié)構(gòu)海明距離，也就是相對于真實的CPDAG錯誤邊數(shù)，包括多出的邊數(shù)(ee)、缺少的邊數(shù)(me)、反向的邊數(shù)(rd)、少方向的邊數(shù)(md,也就是真實情況為有向邊，但算法得到的CPDAG中為無向邊數(shù))、多方向的邊數(shù)(ed,也就是真實情況為無向邊，但算法得到的CPDAG中為有向邊)。

對頂點數(shù)p=20, 樣本量n=100的情況進(jìn)行100次模擬，并利用箱線圖和柱狀圖比較各種算法估計出的模型與真模型之間的結(jié)構(gòu)海明距離，以及各類錯誤邊數(shù)的平均值。模擬結(jié)果分別如圖1和圖2所示。

圖1給出24種算法的SHD箱線圖。可以看到，在BPKL中無論采用GES算法還是PC算法，得到CPDAG的結(jié)構(gòu)海明距離均明顯小于傳統(tǒng)GES和PC算法的結(jié)構(gòu)海明距離。當(dāng)閾值f較小時，與真實圖的結(jié)構(gòu)海明距離也較小，例如gesb10、gesbR10和pcb10、pcbR10。而隨機(jī)的DAG初值，對閾值f的選取依賴更小一些。

圖2給出24種算法的柱狀圖，其中橫軸表示24種算法，縱軸rd-me-md-ee-ed代表在100次模擬中五種錯誤邊出現(xiàn)次數(shù)的均值。在CPDAG的比較結(jié)果中，GDS、GDSM、GES算法的ee最多，PC算法的me最多。在BPKL中，由于在c)中限制了骨架的搜索范圍必須出現(xiàn)在B個DAG中，并且設(shè)置了候選邊的閾值f，因而無論是選用ges，還是選用PC，me的錯誤都增大了一些，特別是對于gesb90，要求候選邊必須出現(xiàn)超過90%的頻率，導(dǎo)致me相對較大。但是這樣的策略，可以大大降低ee、rd、ed錯誤，從而使整體結(jié)構(gòu)海明距離變小。整體來說，gesb10和gesbR10表現(xiàn)最好。

其他條件下的模擬結(jié)果與其類似，因篇幅有限，文中不再給出。需要指出的是,隨著頂點數(shù)p的增加，SHD的整體水平會增加，隨著樣本量n的增加，SHD的整體水平會略微降低。此外，還模擬了B=50,100的情況。我們發(fā)現(xiàn)，B=50或100時,BPKL的表現(xiàn)并不比B=200時差，很難說明B取哪個最優(yōu)。

4 實例分析

利用高斯貝葉斯網(wǎng)建立擬南芥(Arabidopsis thaliana)基因的關(guān)系，該數(shù)據(jù)含有39個基因的118個基因芯片(Affymetrix)微陣列數(shù)據(jù)。Wille等[13]曾利用高斯圖模型刻畫基因間的關(guān)系。先對變量進(jìn)行正態(tài)變換和標(biāo)準(zhǔn)化，以消除各變量之間不同量綱對結(jié)果的影響。然后用模擬分析中的算法對處理后的數(shù)據(jù)進(jìn)行建模。

文中僅展示gesb10、gesbR10、pcb1和pcbR1算法對真實數(shù)據(jù)的估計結(jié)果。通過這4種算法得到的鄰接矩陣(見圖1)可以看到，gesb10和gesbR10的估計結(jié)果相近，pcb1和pcbR1的結(jié)果相差很小，而pcb1和pcbR1的估計結(jié)果明顯稀疏。圖2給出了這4種算法同時估計出來的邊所構(gòu)成的鄰接矩陣，其中橫坐標(biāo)和縱坐標(biāo)代表著每種基因的名字。為使估計結(jié)果更具有可信度，選擇圖2中鄰接矩陣所對應(yīng)的模型作為最終的估計結(jié)果，如圖3所示。

共同的鄰接矩陣如圖4所示。

最終估計圖模型如圖5所示。

5 結(jié) 語

給出了貝葉斯網(wǎng)結(jié)構(gòu)學(xué)習(xí)的新算法BPKL。該方法基于自助法，并利用KL距離衡量最優(yōu)模型到B個DAG的平均距離。通過對不同情況的模擬，發(fā)現(xiàn)自助法樣本量B=50,100,200時，BPKL都有良好的表現(xiàn)。當(dāng)候選邊的閾值f較小時，雖然計算時間稍微增加，但結(jié)構(gòu)學(xué)習(xí)的效果會好很多，可以有效地降低SHD。文中只討論了高斯貝葉斯網(wǎng)的結(jié)構(gòu)學(xué)習(xí)問題,對于離散變量的貝葉斯網(wǎng)和混變量的貝葉斯網(wǎng)依然適用,而且已經(jīng)進(jìn)行了初步嘗試，取得了不錯的模擬結(jié)果[14]。