基于非穩態加性噪聲模型的因果發現算法

郝志峰,丁凱培,蔡瑞初＊,陳薇

(1. 廣東工業大學計算機學院,廣東 廣州 510006;2. 汕頭大學理學院,廣東 汕頭 515063)

0 引言

挖掘可觀測的時間序列之間的因果關系可以幫助人們理解復雜事物的本質和規律。近年來,因果關系被廣泛應用于深度學習[1]、生物信息學[2]、金融經濟[3]、社會科學[4-5]等領域。通常來說,隨機實驗是挖掘因果關系的黃金標準[6],但受限于經濟成本和道德倫理方面的約束,隨機實驗并不總是可行的。因此,如何從可觀測的時間序列中挖掘因果關系成為數據科學的主要問題。

當前,研究人員已經提出了多種恢復時序數據集的因果結構的方法,常見的有基于約束的方法[7-9]、基于因果函數模型的方法[10-11]等。然而,大多數方法都未考慮數據集中非穩態擾動帶來的影響,進而導致方法失效。非穩態擾動引起的數據分布偏移常見于真實的應用環境中,由非穩態擾動所導致的統計誤差會一定程度上打亂數據分布,使數據分布難以表征。一些學者嘗試用數學工具刻畫非穩態擾動的規律,如文獻[12]提出基于約束的異構/非平穩因果發現(CD-NOD)算法,該算法通過引入代理變量從而對非穩態擾動進行表征,并將包含代理變量的增廣變量集作為因果發現算法的輸入。CD-NOD算法首先利用PC(Peter-Clark)算法檢測局部因果機制發生變化的變量,并恢復變量的因果骨架圖,然后通過多項定向規則對因果骨架圖中的無向邊進行定向。然而,CD-NOD算法額外假設數據集中非穩態擾動是彼此獨立的,該假設致使CD-NOD算法無法從不滿足非穩態擾動獨立的數據集中恢復因果結構,并且存在馬爾可夫等價類問題。

針對以上問題,本文基于非穩態擾動與時序信息存在高度相關性的規律,在加性噪聲模型的基礎上將非穩態擾動刻畫成一項關于時序信息的映射并融入模型中,從而提出一種非穩態加性噪聲模型(Nons-ANM)。給出該函數模型的識別條件,并基于模型的識別條件提出一種兩階段的因果發現算法。第1階段基于因果結構圖中葉子節點的噪聲與其余節點均獨立的性質找到因果次序;基于前階段的因果次序先驗,第2階段進一步移除冗余的父子關系,最終得到變量集的因果結構。

1 相關工作

針對觀察數據特性的不同,因果發現方法又可以分為基于時序觀察數據的方法和基于非時序觀察數據的方法[13]。

在非時序方法中,基于約束的方法通常利用條件獨立和忠實性假設來恢復變量之間的因果結構,經典的方法包括文獻[14]介紹的PC算法和文獻[15]介紹的IC(Inductive Causation)算法。這類方法無法解決馬爾可夫等價類的問題,因此算法最終得到的因果結構圖中可能包含未確定方向的邊。為解決馬爾可夫等價類的問題,學者們提出了基于函數的方法。這種方法的核心思想在于:假設存在一種表達力較強的函數形式能夠刻畫變量之間的關系。在這種函數關系下,結果變量能夠由原因變量和噪聲變量的某種映射所表征。在一定的約束條件下,變量之間的因果關系能夠通過檢驗原因變量和噪聲變量的獨立性得出。經典的基于函數的方法包括文獻[16]介紹的非線性加性噪聲模型、文獻[17]介紹的線性非高斯無環模型和文獻[18]介紹的后非線性因果模型等。

時序觀察數據在時間維度上伴隨著“原因在前,結果在后”的特性[13],基于這項特性,非時序方法能夠向時序數據的場景拓展。典型的有文獻[7]提出的PC瞬時條件獨立性(PCMCI)算法,該算法在PC算法的基礎上拓展出恢復多元時序數據集的因果結構的能力。具體來說,PCMCI算法首先使用 PC算法學習得到數據集的馬爾可夫等價類集合,然后再利用瞬時條件獨立性(MCI)檢驗以降低因果關系的誤發現率。然而,PCMCI算法存在未考慮即時效應和隱變量的局限性。因此,文獻[8]提出PCMCI+算法以解決即時效應的問題;文獻[9]提出的LPCMCI算法進一步解決了隱變量的問題。類似地,文獻[10]提出的TiMINo算法和文獻[11]提出的VarLiNGAM算法則是基于函數模型的方法在時序場景上的拓展。

當前針對非穩態數據的因果發現算法的研究仍處于初始階段,文獻[12]介紹的CD-NOD框架著力于在變量的分布發生變化的情況下恢復因果結構,該框架的關鍵部分為引入代理變量C充當潛在變化因素的表征。CD-NOD能夠恢復非穩態數據的因果骨架,但是非穩態擾動的獨立變化假設和馬爾可夫等價類難題限制了該算法對完整因果結構的學習能力。具體來說,馬爾可夫等價類是指具有相同的因果骨架和條件概率分布,但無法確定方向的結構。在CD-NOD算法的定向規則失效后,該算法最終僅能輸出由PC算法學習到的部分定向的因果結構。

2 問題定義

本文所考慮的問題為:給定非穩態時間序列數據集D,如何有效地根據非穩態擾動的特征建模,并恢復非穩態數據集D中n個變量的因果結構。

3 非穩態時序數據的因果結構學習

針對所研究的問題,本節提出一種非穩態加性噪聲模型,并給出該模型在非穩態場景下的可識別性條件及其證明,隨后基于非穩態加性噪聲模型的可識別性理論提出一種兩階段的因果關系發現算法。

3.1 非穩態加性噪聲模型

在穩態場景中,文獻[16]基于二元變量的因果關系提出非線性加性噪聲模型,并證明該模型在限定的函數和噪聲分布下具有可識別性。文獻[10-19]則進一步將非線性加性噪聲模型推廣至時序數據場景和多元變量場景。然而,目前對于非穩態加性噪聲模型的建模及其識別性仍然未知,本文首先給出非穩態加性噪聲模型的定義,隨后給出模型的識別性條件及證明。

定義1(非穩態加性噪聲模型) 給定一組觀測變量集X={X1,X2,…,Xn},若這組觀測變量滿足以下約束,則稱為非穩態加性噪聲模型:

1)觀察變量具有一定的因果次序,因果次序在后的變量不能是次序在前的變量的原因;數據產生過程是遞歸的,變量間的因果關系可以表示為一個有向無環圖。

2)偽因果充分性[12]假設:若存在未觀測到的混淆變量,則這個混淆變量可以用關于時序信息T的光滑函數表示。在任意時刻T=t,這項混淆變量的函數值是確定的。

3)觀測變量Xi∈X服從以下函數關系:

(1)

受文獻[19]提出的多元加性噪聲模型的可識別性條件啟發,Nons-ANM的可識別性理論如定理1所示。在給出定理1之前,先引入(B,F)-IFMOC的識別理論。

引理1(B,F)-IFMOC的識別理論[19-20]:假設一組變量X的聯合分布p(X)是由以因果圖G表征的(B,F)-IFMOC的函數模型導出的,那么這組變量一定不存在兩個不同的(B,F)-IFMOC函數模型表征,同時也一定不能由兩個不同的因果圖G和G′表征。

定理1假設變量集X={X1,X2,…,Xn}滿足非穩態加性噪聲模型,模型屬于(B,F)-IFMOC空間,且增廣變量集X∪{T}的因果結構圖服從因果最小性假設,那么變量集X的因果結構能夠被唯一確定。

證明假設增廣變量集X∪{T}存在兩種不同的非穩態加性噪聲模型表征,分別對應兩個不同的因果結構圖G和G0。

定理1提供了非穩態加性噪聲模型的識別性條件,而對于模型中存在直接因果關系的變量對之間的非對稱性質,可由定理1推論而來。

證明在本文模型中,若兩個變量之間存在因果關系,如Xi→Xj,可分為兩種情況:1)存在時滯的因果關系;2)即時因果關系。兩種情況不會同時發生,因此分別證明在以上兩種情況下本文模型的非對稱性質均成立。

推論1提供了本文模型中原因變量和結果變量之間的非對稱性質。基于上述理論,第3.2節將提出面向非穩態數據的因果發現算法。

3.2 非穩態數據的因果發現算法

基于定理1提供的非穩態加性噪聲模型的識別性條件,本節提出一種兩階段的因果發現算法。算法的初始輸入為觀測變量的時序數據集,最終輸出為表征觀測變量因果結構的有向無環圖,算法框架如圖1所示。

圖1 Nons-ANM算法框架Fig.1 Nons-ANM algorithm framework

算法第1階段主要基于有向無環圖的匯點的性質:對于一個有向無環圖,必然存在出度為0的葉子節點,本文將其稱為“匯點”。在因果結構圖中,若Xsink為匯點,則其噪聲變量為Nsink,它將擁有以下性質:

Nsink⊥{Xk|Xk∈X{Xsink}}

(2)

對于因果結構圖中的匯點Xsink,其噪聲Nsink將獨立于變量集中其余所有變量,這是因為其余的變量Xk均為匯點Xsink的非后裔節點,因此可由不包含匯點Xsink成分的獨立同分布噪聲變量集N〗{Nsink}的某種映射所表征。

基于上述性質,變量的因果次序可通過迭代尋找匯點的方法來獲得。具體來說,對于未確定因果次序的變量Xs∈X,利用變量集中其余未確定因果次序的變量{Xk|Xk∈X{Xs}}對Xs進行回歸并計算得到殘差rs,隨后計算該殘差與其余變量的獨立性得分IND(rs,{Xk|Xk∈X{Xs}}),該得分衡量了變量間的獨立性強弱,衡量方法IND(·)可根據數據性質靈活選擇,本節列舉了一些參考方法。重復上述步驟,選出獨立性得分最高的變量作為匯點,將匯點插入因果次序中,并將匯點的數據從數據集中移除。算法偽代碼如算法1所示。

算法1因果次序學習

輸入觀測變量的時序數據集D={x1,x2,…,xn},最大時延p

輸出觀測變量的因果次序Π

1.初始化變量索引集S={1,2,…,n},因果次序Π為空集

2.WHILE !is Empty(S):

3.FOR s in S:

5.計算IND(rs,{xk|k∈S{s})并記錄

6.END FOR

7.選出獨立性得分最高的變量XM,從變量索引集S中移除M,在因果次序Π中添加M

8.END WHILE

9.RETURN Π

算法1可以學習得到變量集X的因果次序Π,基于因果次序Π可以生成變量集X的部分因果結構。具體來說,在因果次序中,次序在后的變量不能是次序在前的變量的原因。這就意味著,因果次序中包含的確定信息為不存在由次序在后的節點指向次序在前的節點的有向邊;而未確定的信息為次序在前的節點可能存在指向次序在后的節點的有向邊,同時也可能不存在有向邊。因此,因果次序僅表征了因果結構的部分信息,需要進一步消除冗余信息,從而得到最終的因果結構圖。

消除冗余因果關系的步驟如下:

1)對于各節點Xi∈X,以Xi為結果變量,以因果次序在Xi之前的全部節點為原因變量,連接原因變量和結果變量,初始化變量集X的部分因果結構G*。

2)對于G*中可能存在的冗余邊,采用回歸計算得到殘差然后進行獨立性檢驗的方式消除,即對于各節點Xi∈X,依次假設因果次序在Xi之前的某個節點Xj∈{Xj|Xj∈X{Xi},Π(j)<Π(i)}不是Xi的直接原因;在特征數據集中移除Xj的數據后再對Xi的回歸計算得到殘差ri,隨后檢驗殘差ri與變量集{Xj|Xj∈X{Xi},Π(j)<Π(i)}的獨立性。若獨立性仍然成立,則認為Xj不是Xi的直接原因,并移除冗余邊Xj→Xi;否則認為Xj是Xi的直接原因,并保留連接邊。算法偽代碼如算法2所示。

算法2移除冗余因果邊

輸入因果次序Π,觀測變量的時序數據集D={x1,x2,…,xn},最大時延p,顯著性水平α

輸出因果結構圖G

1.基于因果次序Π初始化有向無環圖G*。

2.FOR i in Π:

3.FOR j in {j|Π(j)<Π(i)}:

5.IF IND(ri,{xk|Π(k)<Π(i))>α:

6.刪除G*中的有向邊Xj→Xi,更新父親節點集PAi=PAi{j}

7.END IF

8.END FOR

9.END FOR

10.RETURN G*

根據數據集特征,算法中的回歸方法和獨立性檢驗方法可以相應調整。在因果發現中常見的回歸方法有高斯過程回歸(GPR)、向量自回歸模型(VAR)等;常見的獨立性檢驗方法有偏相關檢驗[22]、希爾伯特-施密特獨立性準則(HSIC)[23-24]、互信息(Ml)[25]、核條件獨立性檢驗(KCI-test)[26]等。本文算法的實際表現受回歸工具和獨立性檢驗工具的影響,統計學的難題會影響算法的實際表現,例如高維變量集、小樣本量等場景。

對本文算法進行復雜度分析:算法1和算法2的時間復雜度均為O(n2·Or(n,l)·Oi(n,l)),其中,n為觀測變量的個數,l為時序數據集的樣本量,Or(n,l)表示算法選用的回歸方法的時間復雜度,Oi(n,l)表示算法選用的獨立性檢驗方法的時間復雜度。

4 實驗

本節將分別基于仿真數據集和真實因果結構數據集開展隨機實驗,從而對本文算法的有效性進行評估。另外采用CD-NOD算法、LPCMCI算法和TiMINo算法作為對比方法,其中CD-NOD算法為非穩態因果關系學習算法,LPCMCI算法是學習含有即時因果效應和隱變量的因果關系學習算法,TiMINo算法則是經典的基于函數因果模型的時序因果關系學習算法。

仿真數據集實驗和真實因果結構數據集實驗中的對比算法均使用官方實現的開源代碼,算法參數均設置為默認值。本文算法Nons-ANM采用高斯過程回歸作為回歸方法,采用基于希爾伯特-施密特獨立性準則的假設檢驗作為獨立性檢驗方法。實驗中所有方法的顯著性水平α均設置為0.05。

在實驗部署環境方面,仿真數據集實驗和真實因果結構數據集實驗均執行于內存為8 GB、CPU型號為Intel i7-7700、操作系統為Windows的環境。

4.1 仿真數據集實驗

在仿真數據集實驗中,本文基于Erdos-Renyi模型隨機生成不同的有向無環圖來表示因果結構,并設計了5組控制變量實驗,具體設置為:變量數量為{5,10, 15, 20},平均入度為{1.0,1.5, 2.0, 2.5},樣本數量為{1 000,2000, 3 000, 4 000},非穩態變量的占比為{0.0,0.4, 0.8, 1.0},最大時延為{0, 2, 4},花括號中的加粗參數為默認的參數設置。變量集中的變量分為穩態和非穩態兩類,對于穩態的變量,其時序數據基于函數因果模型生成,如式(3)所示:

(3)

對于非穩態擾動的變量,其時序數據的生成則基于函數因果模型,如式(4)所示:

(4)

實驗的評價指標選用準確率(P)、召回率(R)和F1值(F1),各評價指標的計算公式如下:

(5)

(6)

(7)

其中:TP表示學習正確的連接邊的數量;FN表示原結構不存在邊,且算法結果中也沒有邊的數量;FP表示原結構不存在邊,但算法結果存在邊的數量。在這套評價算法中,正確率是指學習的因果結構圖中正確的邊數占比;召回率表示正確識別的邊數在所有真實存在的正確邊數中的占比;而F1值結合正確率和召回率,能夠綜合地衡量一種算法的有效性。

本文算法與CD-NOD、LPCMCI和TiMINo算法的實驗結果如圖2所示。

圖2 各算法在不同參數下的實驗結果Fig.2 Experimental results of each algorithms in different parameters

1)圖2(a)～圖2(c)給出不同變量數量下的實驗結果。可以看出,本文算法在所有變量數量的設定下均取得最高的準確率、召回率和F1值,這驗證了本文算法的魯棒性。CD-NOD算法的實驗結果呈現出準確率高于召回率的特點,原因是CD-NOD算法的定向規則中的獨立擾動假設造成了定向誤差。LPCMCI算法擁有從包含隱變量的數據集中恢復因果結構的能力,這項能力幫助LPCMCI算法在實驗默認設置下取得了0.56的F1值,但其總體實驗結果略低于CD-NOD算法,這是因為LPCMCI算法針對隱變量的能力并不完全適用于非穩態擾動的場景。TiMINo算法采用了保守的策略,當算法無法確定因果關系時將執行熔斷機制;由圖3(a)可知,TiMINo算法在實驗默認設置下熔斷率為0.8;而當變量數量增至20時,其熔斷率達到1.0,即TiMINo算法無法恢復任何一組仿真數據集的因果結構。由此可知,經典的時序因果發現方法難以直接應用于非穩態場景。

圖3 TiMINo算法在不同參數設置下的未識別率Fig.3 Unrecognition rate of TiMINo algorithm under different parameter settings

2)從圖2(d)～圖2(f)可以看出,實驗中的所有算法對樣本數量都比較敏感。在樣本數量逐漸增加的過程中,本文算法的F1值逐漸提升,這是因為高樣本數量能降低統計誤差,從而幫助本文算法逼近理論效果。在所有樣本數量的實驗設置下,本文算法均取得了最高的F1值,這體現了本文算法的優越性。

3)圖2(g)～圖2(i)展示了不同平均入度下的實驗結果。可以看出,CD-NOD算法和LPCMCI算法對平均入度較敏感,兩個算法的F1值隨著平均入度增加而明顯下降,這表明稠密的因果結構更易于暴露算法的不足。TiMINo算法的熔斷率在所有設定中均不小于0.7,這表明了TiMINo算法對非穩態數據集的學習能力不足。相比之下,本文算法在所有設置下的均表現穩定且優勢明顯,這表示本文算法能夠有效地學習稠密的因果結構。

4)圖2(j)～圖2(l)給出了不同非穩態節點占比的實驗結果。可以看出,在所有算法中本文算法擁有最好的非穩態因果關系學習能力,即使在非穩態節點占比為1.0時,本文算法依然能取得0.81的F1值。隨著非穩態擾動節點占比的增加,其余3個算法的實驗結果均大幅降低,這進一步體現了本文算法的魯棒性和優越性。

5)圖2(m)～圖2(o)給出了不同最大時延下的實驗結果。最大時延的增加會增大時序因果網絡的復雜性,進而增大因果發現算法的學習難度。從圖中可以看出,CD-NOD算法、LPCMCI算法和TiMINo算法在最大時延為4時難以取得有效的實驗結果,這是因為以上3種算法均沒有學習復雜的時序因果網絡的能力。而本文算法在所有最大時延的實驗設置下表現穩定,實驗結果顯著優于其他3種算法,說明了本文算法的有效性。

6)圖3表示TiMINo算法在不同控制實驗中觸發熔斷策略的數據集比率,即TiMINo算法的未識別率。可以看出,TiMINo算法頻繁觸發了熔斷機制:在所有控制實驗中,TiMINo的未識別率均不小于0.7;而在部分復雜的因果結構中的未識別率達到1.0,例如20個節點和非穩態節點占比達到1.0等。這進一步表明,未考慮非穩態擾動的時序因果發現算法難以直接應用于非穩態場景中。

4.2 真實因果結構數據集實驗

本節采用3種不同規格的真實因果結構進行實驗,相關的真實因果結構信息可參考開源網站:https:∥www.bnlearn.com/bnrepository/。本節實驗中選用的真實因果網絡的參數表1如所示。

表1 真實因果網絡結構參數Table 1 Real causal network structure parameters

本節實驗中采用的算法參數和實驗參數均與仿真數據集實驗相同,同樣采用準確率、召回率以及F1值作為算法效果的評價指標。

各算法的實驗結果如圖4所示。由圖4可知,本文算法在所有數據集中均取得了最優的實驗結果。特別地,本文算法在sachs結構的數據集中取得的F1值為0.86,相比CD-NOD的0.52提升了65.4%,相比LPCMCI的0.53提升了62.3%,這是因為sachs結構比較稠密,而本文算法在稠密的因果結構下更加能夠體現其優越性。而對于變量個數少且相對稀疏的因果結構(如survey結構),本文算法和其他對比算法的F1值接近,這則是因為簡單的因果結構難以體現非穩態擾動造成的影響。總體而言,真實因果結構數據集的實驗結果驗證了本文算法的有效性。TiMINo算法在不同真實結構下的未識別率如圖5所示。

圖4 不同算法的真實結構實驗結果Fig.4 Experimental results of real structures of different algorithms

圖5 TiMINo算法在不同真實結構下的未識別率Fig.5 Unrecognition rate of TiMINo algorithm under different real structures

5 結束語

現有因果發現方法難以恢復非穩態數據集的因果結構,本文提出非穩態加性噪聲模型,并基于該模型的識別性理論提出一種魯棒的兩階段算法,從而為非穩態數據集的因果結構學習難題提供解決方案。為驗證本文算法有效性,引入CD-NOD算法、LPCMCI算法和TiMINo算法,分別在仿真數據集和真實結構數據集上部署實驗,最終的實驗結果體現了本文算法相較于主流因果發現算法在準確率和魯棒性等方面的優勢。下一步將對模型的加性噪聲假設進行弱化并優化算法流程,以擴展本文方法的適用范圍。