李洪飛，萬亞平,2，陽小華,2，耿家興

(1.南華大學(xué) 計(jì)算機(jī)學(xué)院，湖南 衡陽 421001；2.中核集團(tuán)高可信計(jì)算重點(diǎn)學(xué)科實(shí)驗(yàn)室，湖南 衡陽 421001)

0 引 言

因果發(fā)現(xiàn)旨在開發(fā)從觀測中學(xué)習(xí)因果網(wǎng)絡(luò)結(jié)構(gòu)的算法，而因果結(jié)構(gòu)學(xué)習(xí)是機(jī)器學(xué)習(xí)和統(tǒng)計(jì)領(lǐng)域的新課題之一[1-2]。包括預(yù)防科學(xué)在內(nèi)的許多實(shí)證科學(xué)中，需要研究各種現(xiàn)象背后的因果機(jī)制。

通常，從觀測數(shù)據(jù)中發(fā)現(xiàn)因果結(jié)構(gòu)是困難的，因?yàn)榭赡艿哪Ｐ偷某笖?shù)集，造成一些模型可能無法與其他模型區(qū)分開來。在變量p數(shù)量相當(dāng)甚至超過樣本n數(shù)量的情況下是復(fù)雜的。盡管該問題本身存在困難，但已經(jīng)開發(fā)了許多算法。目前推斷算法一般歸為兩類：貝葉斯結(jié)構(gòu)學(xué)習(xí)算法[3-5]和基于加性噪聲模型因果發(fā)現(xiàn)算法[6-9]。例如PC算法[3]是一種基于約束的貝葉斯結(jié)構(gòu)學(xué)習(xí)算法，它先推斷出一組條件獨(dú)立關(guān)系，然后識別相關(guān)的馬爾可夫等價(jià)類。但這類方法無法區(qū)分X→Z→Y和X←Z←Y兩種結(jié)構(gòu)。加性噪聲模型中，在各種假設(shè)下，可從觀測數(shù)據(jù)中恢復(fù)準(zhǔn)確的圖像[10]。如Shimizu等假設(shè)數(shù)據(jù)遵循線性非高斯無環(huán)模型(LiNGAM)[7]，后續(xù)的DirectLiNGAM[8]和PLiNGAM[11]方法通過兩兩統(tǒng)計(jì)量迭代選擇因果順序。LiNGAM提出的這些方法在高維設(shè)置中是不一致的，這些設(shè)置允許變量數(shù)的伸縮與樣本大小相同或更快，并且只在線性模型下適用。Hoyer針對非線性數(shù)據(jù)因果模型，提出適用于連續(xù)數(shù)據(jù)算法ANM[6]。為擴(kuò)展到離散數(shù)據(jù)的情況，對ANM算法進(jìn)一步改進(jìn)，提出后非線性數(shù)據(jù)的PostNonLinear算法[9]和信息幾何方法[12]。

近年來，Makoto等提出最小二乘獨(dú)立回歸算法(LSIR)[13]。該算法通過交叉驗(yàn)證自然地優(yōu)化內(nèi)核寬度和正則化參數(shù)等調(diào)優(yōu)參數(shù)，從而避免以依賴數(shù)據(jù)的方式過度擬合。這些方法都只考慮了少量的變量，一旦在高維的情況下(n>7)，因果發(fā)現(xiàn)方法的輸出可能高度依賴于順序，準(zhǔn)確率會(huì)很低。在大多數(shù)情況下是兩個(gè)變量進(jìn)行方向識別。因?yàn)楦鶕?jù)理論和算法，可擴(kuò)展性有限，所以這些因果推斷方法都不適用在高維數(shù)據(jù)的網(wǎng)絡(luò)結(jié)構(gòu)中。

為解決在高維數(shù)據(jù)下學(xué)習(xí)到較準(zhǔn)確的因果網(wǎng)絡(luò)結(jié)構(gòu)問題，曾千千等提出基于MIC和MI的因果推斷算法[14-15]。這些算法將高維網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)問題分解成網(wǎng)絡(luò)中每個(gè)節(jié)點(diǎn)對應(yīng)的低維因果網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)問題。這些算法中，基于信息論的互信息(MI)或者最大信息系數(shù)(MIC)能夠較好地刪除不相關(guān)的節(jié)點(diǎn)，降低計(jì)算復(fù)雜度，但MIC對于高維變量不可用，對于大樣本問題耗時(shí)MI也不穩(wěn)定。隨著樣本量的增大，方差沒有變小，對異常值的魯棒性較弱。由Jiang等提出的CDC相比較互信息和最大信息系數(shù)，可更為準(zhǔn)確地檢測出變量之間的關(guān)聯(lián)關(guān)系[16]。

基于此，文中提出一種基于CDC的高維數(shù)據(jù)因果推斷算法。該算法利用CDC對變量間的依賴度進(jìn)行檢測，再利用條件獨(dú)立性檢測精煉目標(biāo)節(jié)點(diǎn)的父子節(jié)點(diǎn)集，然后使用非線性最小二乘獨(dú)立回歸(LSIR)(兩個(gè)變量之間效應(yīng)方向的新測度)，為圖中的目標(biāo)節(jié)點(diǎn)及其父子節(jié)點(diǎn)之間的無向邊標(biāo)注因果方向，最后迭代所有的節(jié)點(diǎn)完成完整的因果網(wǎng)絡(luò)結(jié)構(gòu)。由于CDC對大樣本下異常值的魯棒性強(qiáng)，因此提高了算法的準(zhǔn)確率。實(shí)驗(yàn)也表明在高維數(shù)據(jù)下該算法的精確度優(yōu)于其他因果推斷算法。

1 預(yù)備知識

1.1 加性噪聲模型

1.2 最小二乘獨(dú)立回歸

最小二乘獨(dú)立回歸(LSIR)是通過最小化輸入和殘差之間的平方損失互信息估計(jì)值來學(xué)習(xí)加性噪聲模型，識別兩個(gè)非線性變量的因果關(guān)系的方法。LSIR相對現(xiàn)有方法的一個(gè)顯著優(yōu)勢是，調(diào)優(yōu)參數(shù)(如內(nèi)核寬度和正則化參數(shù))可以通過交叉驗(yàn)證自然地進(jìn)行優(yōu)化，從而避免以依賴數(shù)據(jù)的方式過度擬合。LSIR與最先進(jìn)的因果推理方法相比是比較有利的。

(1)

采用平方損失相互信息(SMI)[14]作為獨(dú)立措施:

(2)

使用SMI，獲得一個(gè)分析形式的估計(jì)值。得到SMI估計(jì)量后，開始學(xué)習(xí)參數(shù)β回歸模型：

這種方法稱為最小二乘獨(dú)立回歸。

對于回歸參數(shù)學(xué)習(xí)，可以簡單地采用梯度下降法:

(4)

其中，η是步長，可以選一些近似線搜索方法，如Armijo規(guī)則[14]。

(5)

(6)

為了確定因果關(guān)系的方向，計(jì)算兩個(gè)方向X→Y(即X引起Y)和X←Y(即Y導(dǎo)致X)的p值pX→Y和pX←Y。對于一個(gè)給定的顯著性水平δ，確定的因果方向如下。

如果pX→Y>δ和pX←Y≤δ，模型選擇X→Y；

如果pX←Y>δ和pX←Y≤δ，X←Y模型被選中；

如果pX→Y,pX←Y≤δ，X和Y之間沒有因果關(guān)系；

如果pX→Y,pX←Y>δ,則建模假設(shè)是不正確的。

總之，檢驗(yàn)因果模型Y=fY(X)+EY或替代X=fX(Y)+EX是否與數(shù)據(jù)吻合較好，其中擬合優(yōu)度是通過輸入與殘差之間的獨(dú)立性來衡量的(即噪聲估計(jì))，輸入和殘差的獨(dú)立性可以通過置換測試在實(shí)踐中確定[14]。

1.3 Copula相關(guān)系數(shù)

Copula相關(guān)系數(shù)(CDC)用來測量兩個(gè)隨機(jī)向量X和Y之間的依賴關(guān)系，是一種較好的關(guān)聯(lián)檢測魯棒依賴測度。CDC對異常值更具有魯棒性，適用于更廣泛的模型，尤其適用于高維問題。

定理1(概率積分變換)：對于具有累積分布函數(shù)(CDF)F的連續(xù)隨機(jī)變量X，隨機(jī)變量U=F(X)在[0,1]上均勻分布。因此，向量：

U=(U1,U2,…,Up)=P(X)=(F1(X1),F2(X2),…,Fp(Xp))

(7)

稱為Copula變換，具有均勻的邊緣。

定理2：讓隨機(jī)向量X=(X1,X2,…,Xp)連續(xù)邊際CDFs,Fi,1≤i≤p。然后X的聯(lián)合累積分布函數(shù)F(X),X=(x1,x2,…,xp)，唯一表示為:

F(X)=Cx(F1(x1),F2(x2),…,Fk(xp))

(8)

其中分布函數(shù)Cx稱為X的Copula。

定理3：設(shè)X和Y為連續(xù)隨機(jī)變量，那么X和Y是獨(dú)立的，當(dāng)且僅當(dāng)CXY(FX(x),FY(y))=FX(x)FY(y)，其中FX(x)和FY(y)分別為x和y的分布函數(shù)。

首先考慮特殊情況p=q=1。如果X和Y是獨(dú)立的，將FX(X)andFY(Y)分別表示為X和Y的分布函數(shù)，由定理3可知：

CXY(FX(x),FY(y))=FX(x)FY(y)

(9)

如果X和Y是獨(dú)立的，有W=FX(X)andV=FY(Y)是獨(dú)立的。因此，將多維依賴檢驗(yàn)問題轉(zhuǎn)化為零假設(shè)(H0)，W和V是獨(dú)立的二維依賴檢驗(yàn)問題，可以通過MIC、互信息等提出的另一種依賴測度來求解。

定義1：設(shè)X和Y為兩個(gè)隨機(jī)變量，X和Y的最大相關(guān)系數(shù)(MCC)為：

(10)

如果將φ1,φ2限制為線性函數(shù)，MCC是皮爾森相關(guān)系數(shù)。為了計(jì)算MCC，通常會(huì)將限制φ1,φ2限制為再生核希爾伯特空間理論(RKHS)，例如n階的多項(xiàng)式函數(shù)的空間。而ACE用來計(jì)算MCC。

定義2：利用上述表示法，兩個(gè)隨機(jī)向量X和Y之間的CDC由FX(X)andFY(Y)之間的MCC給出：

CDC(X,Y)=MCC(FX(X),FY(Y))

(11)

其中，F(xiàn)X(x)andFY(y)分別為X和Y的分布函數(shù)。

可以看出，CDC是一個(gè)基于秩的依賴性測度，因?yàn)榉植己瘮?shù)是基于秩的，所以CDC適用于高維變量。在實(shí)際應(yīng)用中，真實(shí)的邊際分布函數(shù)是未知的，用經(jīng)驗(yàn)邊際分布或估計(jì)邊際分布代替，得到CDC的估計(jì)。

2 因果結(jié)構(gòu)學(xué)習(xí)算法

本節(jié)根據(jù)相較于互信息和最大信息系數(shù)更為準(zhǔn)確地檢測出變量間關(guān)聯(lián)關(guān)系的CDC，結(jié)合發(fā)現(xiàn)因果骨架的結(jié)構(gòu)學(xué)習(xí)，二元變量的方向識別，最終迭代得到一個(gè)完整表現(xiàn)出高維數(shù)據(jù)間因果關(guān)系的因果網(wǎng)絡(luò)結(jié)構(gòu)。研究表明，張等基于互信息與曾等基于最大信息系數(shù)的算法構(gòu)成的無向圖，對大樣本和高維下異常值的魯棒性不高，一定程度會(huì)增加父子節(jié)點(diǎn)集，造成條件獨(dú)立測試的計(jì)算復(fù)雜度，得到的最終網(wǎng)絡(luò)與真實(shí)網(wǎng)絡(luò)差異很大[14-15]。相較于MIC和MI，CDC能夠更好地度量節(jié)點(diǎn)間的依賴關(guān)系，減少冗余帶來的計(jì)算復(fù)雜度，從而提高效率，并獲得更加準(zhǔn)確的因果網(wǎng)絡(luò)結(jié)構(gòu)。該算法基本框架如圖1所示。

圖1 算法的基本框架

2.1 無向子圖構(gòu)造

CDC可以快速對隨機(jī)變量之間的依賴關(guān)系進(jìn)行檢測。在樣本量增多的情況下，CDC是最佳的依賴測度，它對異常值具有魯棒性，計(jì)算效率高，對大多數(shù)函數(shù)類型具有較好的適用性。對隨機(jī)變量X和Y計(jì)算CDC(X,Y)，X和Y之間的依賴程度與CDC的值成正比。當(dāng)CDC值越大，則兩個(gè)變量對應(yīng)的點(diǎn)相連程度越高。反之，若X和Y之間CDC值很小，甚至達(dá)到0，則兩個(gè)變量之間是獨(dú)立的，對應(yīng)的兩點(diǎn)不相連。根據(jù)CDC的優(yōu)點(diǎn)計(jì)算變量間的CDC(X,Y)值構(gòu)造無向圖，步驟如下：

(1)對數(shù)據(jù)集中任意的兩個(gè)隨機(jī)變量計(jì)算其之間的CDC值，找到每個(gè)節(jié)點(diǎn)與其他節(jié)點(diǎn)的CDC值中最大的值，用MCDC(X,Y)表示；

(2)對于任意節(jié)點(diǎn)，假設(shè)為y，其他節(jié)點(diǎn)表示為X集，如果滿足下列條件，則直接連接在變量X和Y之間建立：CDC(xi,y)≥α·MCDC(X)，否則將xi從X中移除；

(3)將任意的xi∈X加入到PC(y)，應(yīng)用條件獨(dú)立性測試(CI)刪除PC(y)中非父子節(jié)點(diǎn)；

(4)重復(fù)步驟3，迭代完X中所有節(jié)點(diǎn)。

上述步驟中，0.5≤α≤1是一個(gè)參數(shù)，來確定此階段的連接數(shù)量，如果該參數(shù)非常高(接近1.0)，在該算法的這個(gè)階段，存在多個(gè)真邊被拒絕的高概率。另一方面，將該參數(shù)設(shè)置的低，可能會(huì)導(dǎo)致在算法的這個(gè)階段包含幾個(gè)錯(cuò)誤的邊。對于該研究中的問題，參數(shù)α設(shè)置為0.9，并且發(fā)現(xiàn)包括大多數(shù)真實(shí)邊緣，同時(shí)允許在網(wǎng)絡(luò)結(jié)構(gòu)中僅添加少量錯(cuò)誤邊緣。

2.2 目標(biāo)節(jié)點(diǎn)和父子節(jié)點(diǎn)集PC(y)中每個(gè)節(jié)點(diǎn)進(jìn)行方向識別

在前一階段用CDC、CI測試得到的父子節(jié)點(diǎn)集，構(gòu)成了一個(gè)準(zhǔn)確的骨架。再利用LSIR算法對骨架之間每兩個(gè)節(jié)點(diǎn)之間的無向邊進(jìn)行方向識別。相比較于ANM、IGCI二元算法，LSIR在大樣本下通過交叉驗(yàn)證進(jìn)行優(yōu)化，避免了過分依賴數(shù)據(jù)而導(dǎo)致的過度擬合。雖然LSIR算法不能有效處理高維數(shù)據(jù)下的因果網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)問題，但由于在算法2.1無向圖構(gòu)造過程中，已經(jīng)得到了目標(biāo)節(jié)點(diǎn)y的父子節(jié)點(diǎn)集合PC(y)，所以可以利用LSIR對PC(y)中的每一個(gè)變量和y之間進(jìn)行方向判別，這等同于在2維間應(yīng)用LSIR。這種分治策略保證了它的有效性。

2.3 算法描述

文中算法記為CDC & LSIR causal inference (CLCI)。

Input:variable set X={x1,x2,…,xn};threshold αOutput:DAG G1 for each xi∈X doSet xi=y,xi∈X,PC(y)={}/? 選取X中任意一個(gè)節(jié)點(diǎn)為目標(biāo)節(jié)點(diǎn)?/2 for each xi∈X do /? 對PC(y)進(jìn)行精煉簡化,除掉非父節(jié)點(diǎn)?/if CDC(y;xi)<α?MCDC(X), then X=Xxi3 for each xi∈X do /? 構(gòu)造關(guān)于節(jié)點(diǎn)y的無向圖?/ add xi to PC(y),if there is a set S ( arbitrary subset of PC(y)xi),such that y⊥xi|S,then PC(y)=PC(y)xi4 for each xi∈PC(y) doif there is a set S(arbitrary subset of PC(y)xi), such that y⊥xi|S, then PC(y)=PC(y)xi5 for each xi∈PC(y) do /?對骨架之間每兩個(gè)點(diǎn)的方向進(jìn)行判別? /employs LSIR algorithm to distinguish the direction of (y,xi)

3 模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)設(shè)置與環(huán)境

通過人工合成模擬數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)一用于在大樣本下測試CDC的性能，而實(shí)驗(yàn)二用于驗(yàn)證CLCI算法的有效性。實(shí)驗(yàn)中盡可能假設(shè)復(fù)雜條件，即觀察樣本數(shù)目足夠大，并且模糊了對噪聲變量非高斯性的假設(shè)。每個(gè)節(jié)點(diǎn)在因果網(wǎng)絡(luò)結(jié)構(gòu)中的序列按照函數(shù)y=ω1*tan(x1)+ω2*tan(x2)+ε生成。每個(gè)函數(shù)中的權(quán)值分別為ω1、ω2,隨機(jī)取值0.3～0.7之間，y的父節(jié)點(diǎn)分別為X1和X2。噪聲ε服從均勻分布，且權(quán)值為6%。文中提出的算法記為CLCI。實(shí)驗(yàn)平臺為Window 7 64 bit操作系統(tǒng)，配置為Intel酷睿i5-4200U，內(nèi)存4 GB，實(shí)驗(yàn)環(huán)境為Matlab-R2016a。

3.2 實(shí)驗(yàn)分析

在不同樣本下，CDC與MI和MIC計(jì)算時(shí)間對比如圖2所示(表示時(shí)間效率)。

圖2 MI、MIC和CDC在不同大樣本下的計(jì)算時(shí)間對比

從結(jié)果可以看出，隨著樣本量的增大，MIC的運(yùn)行時(shí)間迅速增加，MI直到樣本大小約為4 500時(shí)才開始緩慢增加，而CDC的運(yùn)行時(shí)間則沒有增加趨勢。

為評價(jià)CLCI在高維下的優(yōu)劣，避免隨機(jī)性，在實(shí)驗(yàn)中選擇5 000樣本量分別生成20維到200維的因果網(wǎng)絡(luò)圖，用來測試高維度下不同樣本量的實(shí)驗(yàn)效果。實(shí)驗(yàn)引入三個(gè)常用指標(biāo)recall、precision和F1來評價(jià)算法性能。在不同維度下四種算法的評分參數(shù)如圖3所示。

從圖3可以看到，IGCI算法的召回率在超過40維左右時(shí)比其他三種算法都要高，但由于在高維因果網(wǎng)絡(luò)結(jié)構(gòu)下，IGCI錯(cuò)誤地添加了很多邊，準(zhǔn)確率相對較低。同樣，加性噪聲模型(ANM)算法在高維數(shù)據(jù)情況下三個(gè)評分參數(shù)比CLCI和TPCI算法低。由于CLCI和TPCI算法采用了降維的方法，在維數(shù)增加的過程中，其保持了較好的穩(wěn)定性，但是相對于利用MI的TPCI，CLCI利用CDC在大樣本下的魯棒性更好，時(shí)間復(fù)雜度更低，耗時(shí)更短。由圖3也可以看出，CLCI要優(yōu)于其他三種算法。

(a)recall

(b)precision

(c)F1

4 結(jié)束語

與其他適用于高維數(shù)據(jù)的因果推斷方法不同，文中結(jié)合基于信息論的Copula依賴系數(shù)，可以在高維和大樣本數(shù)據(jù)下更為準(zhǔn)確地檢測出變量之間的關(guān)聯(lián)關(guān)系，降低計(jì)算復(fù)雜度，再結(jié)合條件獨(dú)立性測試對數(shù)據(jù)集進(jìn)行無向結(jié)構(gòu)學(xué)習(xí)，最后用LSIR算法對無向結(jié)構(gòu)骨架進(jìn)行每兩點(diǎn)間的方向判斷，得到最終的因果網(wǎng)絡(luò)結(jié)構(gòu)。文中采用虛擬網(wǎng)絡(luò)進(jìn)行的實(shí)驗(yàn)表明，利用CDC進(jìn)行無向骨架的構(gòu)造，耗時(shí)短，計(jì)算復(fù)雜度低于其他算法。當(dāng)數(shù)據(jù)集維數(shù)較高、樣本量大時(shí)，利用分治策略，該算法要大大優(yōu)于其他因果推斷算法。