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

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

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

0 引 言

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

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

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

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

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

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

b)得分搜索階段，在這些條件獨(dú)立性關(guān)系的約束下，搜索評(píng)分最低的DAG。例如，Tsamardinos等[4]提出MMHC算法，先采用基于約束的方法找到一個(gè)頂點(diǎn)可能的父頂點(diǎn)和子頂點(diǎn)，然后搜索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ù)集計(jì)算穩(wěn)健的模型選擇準(zhǔn)則。對(duì)于每個(gè)數(shù)據(jù)集和一個(gè)給定的DAG，求此DAG參數(shù)的極大似然估計(jì)，進(jìn)而利用聚合的信息準(zhǔn)則來(lái)衡量G;Wang 等[10]提出了一種基于Bootstrap的DAG結(jié)構(gòu)學(xué)習(xí)方法，稱為DAGBag。該方法基于數(shù)據(jù)的Bootstrap重采樣學(xué)習(xí)DAG集合，以結(jié)構(gòu)海明距離(SHD)為度量，通過(guò)最小化到整個(gè)DAG集合的總距離來(lái)導(dǎo)出聚合DAG。

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

1)對(duì)觀測(cè)樣本重抽樣得到B組Bootstrap樣本；

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

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

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

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

貝葉斯網(wǎng)結(jié)構(gòu)由有向無(wú)圈圖G(V,E)表示。頂點(diǎn)集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ù)可以分解為

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

當(dāng)聯(lián)合分布P為多元高斯分布時(shí)，高斯貝葉斯網(wǎng)可根據(jù)密度函數(shù)的分解寫(xiě)成結(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維單位矩陣，通過(guò)簡(jiǎn)單的矩陣變換有

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，

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

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

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

樣本均值

樣本協(xié)方差陣

則有對(duì)數(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òng)，都可能導(dǎo)致得到的DAG有較大變化。為此，Wang 等[10]提出了自助聚合(Bootstrap aggregating)方法，以得到更穩(wěn)定的圖結(jié)構(gòu)。這種方法可以在一定程度上避免過(guò)擬合，并減小學(xué)習(xí)算法的假邊發(fā)現(xiàn)率。它們的算法稱為DAGBag，詳細(xì)過(guò)程見(jiàn)表1。

表1 算法1. DAGBag

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

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

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

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

從而有

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

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

選取KL距離作為度量，是因?yàn)椋?/p>

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

以及正定矩陣Σ，有

其中

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

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

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

在算法3中，有以下幾點(diǎn)需要說(shuō)明。首先，針對(duì)算法c)中的優(yōu)化問(wèn)題

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

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

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

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

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

3 模擬實(shí)驗(yàn)

通過(guò)模擬來(lái)比較BPKL算法、GDS算法、基于道義圖的GDSM算法、PC算法和GES算法。

在BPKL算法中，a)重抽樣次數(shù)設(shè)置為B=200;b)結(jié)構(gòu)學(xué)習(xí)算法L采用GES算法和PC算法;c)搜索過(guò)程中，骨架候選邊頻率的閾值，對(duì)于GES算法，選為f=10%,30%,50%,70%,90%，而對(duì)于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的空間中搜索的，因而得到一個(gè)DAG；而GES和PC將得到CPDAG，用于表示馬爾可夫等價(jià)。為了比較公平性，將前三種方法得到的DAG轉(zhuǎn)換成CPDAG。從而將它們的CPDAG求與真實(shí)模型的CPDAG之間的結(jié)構(gòu)海明距離。所謂結(jié)構(gòu)海明距離，也就是相對(duì)于真實(shí)的CPDAG錯(cuò)誤邊數(shù)，包括多出的邊數(shù)(ee)、缺少的邊數(shù)(me)、反向的邊數(shù)(rd)、少方向的邊數(shù)(md,也就是真實(shí)情況為有向邊，但算法得到的CPDAG中為無(wú)向邊數(shù))、多方向的邊數(shù)(ed,也就是真實(shí)情況為無(wú)向邊，但算法得到的CPDAG中為有向邊)。

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

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

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

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

4 實(shí)例分析

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

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

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

最終估計(jì)圖模型如圖5所示。

5 結(jié) 語(yǔ)

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