蔡瑞初，伍運(yùn)金，陳薇，郝志峰，2

（1.廣東工業(yè)大學(xué) 計(jì)算機(jī)學(xué)院，廣州 510006；2.汕頭大學(xué) 理學(xué)院，廣東 汕頭 515063）

0 概述

因果發(fā)現(xiàn)旨在從觀測(cè)數(shù)據(jù)中發(fā)現(xiàn)變量之間的因果關(guān)系，可以揭露數(shù)據(jù)的生成機(jī)制，幫助人們理解數(shù)據(jù)，輔助人們進(jìn)行干預(yù)和決策［1］。近年來(lái)，因果關(guān)系在深度學(xué)習(xí)［2］、金融經(jīng)濟(jì)［3］、神經(jīng)科學(xué)［4］、生物信息學(xué)［5］、社會(huì)科學(xué)［6］等領(lǐng)域受到了廣泛關(guān)注。

當(dāng)觀測(cè)數(shù)據(jù)是多元時(shí)間序列數(shù)據(jù)時(shí)，現(xiàn)有的時(shí)序因果發(fā)現(xiàn)算法［7-10］通常認(rèn)為個(gè)體之間是獨(dú)立的，為每一個(gè)個(gè)體的多元時(shí)間序列數(shù)據(jù)單獨(dú)學(xué)習(xí)一個(gè)因果關(guān)系作為該個(gè)體背后的因果關(guān)系，而個(gè)體間的因果關(guān)系學(xué)習(xí)過程是彼此無(wú)關(guān)的。然而在實(shí)際中，個(gè)體之間可能存在相同的因果關(guān)系。舉例來(lái)說，在電商場(chǎng)景下，同一群體中的個(gè)體的購(gòu)買行為可能具有相同的思維方式，如對(duì)于家庭群體的個(gè)體是否購(gòu)買某個(gè)商品的影響變量是商品的質(zhì)量，而對(duì)于普通家庭群體的個(gè)體是否購(gòu)買某個(gè)商品的影響變量則是商品的質(zhì)量和價(jià)格。

因此，來(lái)源于不同群體的個(gè)體的數(shù)據(jù)背后會(huì)具有不同的因果關(guān)系（產(chǎn)生機(jī)制），而相同群體的個(gè)體的數(shù)據(jù)背后會(huì)具有相同的因果關(guān)系。如果能知道哪些個(gè)體屬于同一群體，便能利用群體內(nèi)多個(gè)獨(dú)立同分布的個(gè)體數(shù)據(jù)一起學(xué)習(xí)該群體共同的因果關(guān)系。然而，在實(shí)際中無(wú)法預(yù)先知道不同個(gè)體是否屬于相同群體，也無(wú)法判斷總共有多少個(gè)群體。

針對(duì)上述問題，本文提出一種面向多元時(shí)間序列數(shù)據(jù)的群體因果關(guān)系發(fā)現(xiàn)算法。首先基于因果關(guān)系的相似性，將所有個(gè)體劃分成多個(gè)群體且無(wú)須指定群體的個(gè)數(shù)。對(duì)于每一個(gè)群體，使用變分推斷方法學(xué)習(xí)群體因果關(guān)系，從而充分利用多個(gè)個(gè)體數(shù)據(jù)。當(dāng)所有個(gè)體均在一個(gè)群體時(shí)，該算法利用所有個(gè)體數(shù)據(jù)共同學(xué)習(xí)一個(gè)因果關(guān)系。當(dāng)一個(gè)群體內(nèi)只有一個(gè)個(gè)體時(shí)，該算法與現(xiàn)有時(shí)序因果關(guān)系發(fā)現(xiàn)算法類似，僅利用單個(gè)個(gè)體數(shù)據(jù)單獨(dú)學(xué)習(xí)一個(gè)因果關(guān)系。

1 相關(guān)工作

因果關(guān)系發(fā)現(xiàn)算法按照觀測(cè)數(shù)據(jù)的類型可以分為基于非時(shí)序數(shù)據(jù)的因果發(fā)現(xiàn)方法和基于時(shí)序數(shù)據(jù)的因果發(fā)現(xiàn)方法［11］。

基于非時(shí)序數(shù)據(jù)的因果發(fā)現(xiàn)方法中包括基于約束的方法［12］、基于評(píng)分的方法［13］和基于函數(shù)的方法［14-16］?；诩s束的方法利用（條件）獨(dú)立性檢驗(yàn)來(lái)判斷變量之間是否存在因果關(guān)系，而基于評(píng)分的方法通過給DAG 打分并尋找得分最高的DAG 作為變量間的因果關(guān)系，但2 種方法都存在馬爾可夫等價(jià)類的問題。為了解決這個(gè)問題，學(xué)者們提出了基于函數(shù)的方法，此類方法從數(shù)據(jù)產(chǎn)生機(jī)制出發(fā)，假設(shè)原因變量與結(jié)果變量存在函數(shù)映射，以及存在與原因變量獨(dú)立的噪聲變量，通過原因與噪聲的獨(dú)立性來(lái)識(shí)別因果關(guān)系?；诤瘮?shù)的方法包括線性非高斯無(wú)環(huán)模型（Linear Non-Gaussan Acyclic Model，LiNGAM）［14］、非線性加性噪聲模型（nonlinear Additive Noise Model，ANM）［15］和后非線性因果模型（Post-NonLinear causal model，PNL）［16］。

上述基于非時(shí)序數(shù)據(jù)的因果發(fā)現(xiàn)方法也被拓展到了時(shí)序數(shù)據(jù)上，如同樣基于約束的PCMCI算法［8］、基于評(píng)分的DYNOTEARS算法［9］和基于函數(shù)的VAR-LiNGAM算法［10］。PCMCI 算法基于條件獨(dú)立性測(cè)試框架，使用PC 算法［17］發(fā)現(xiàn)變量的馬爾可夫等價(jià)類集合從而縮短條件集，進(jìn)一步用瞬時(shí)條件獨(dú)立性（Momentary Conditional Independence，MCI）檢驗(yàn)降低誤發(fā)現(xiàn)率。DYNOTEARS算法從優(yōu)化問題的角度出發(fā)，最小化一個(gè)帶無(wú)環(huán)約束的損失函數(shù)，從而學(xué)習(xí)變量之間的瞬時(shí)影響和時(shí)延影響。VAR-LiNGAM 算法結(jié)合自回歸模型和線性非高斯無(wú)環(huán)模型來(lái)識(shí)別變量的瞬時(shí)和時(shí)延因果關(guān)系的影響權(quán)重，并通過非高斯性假設(shè)保證了算法的可識(shí)別性。

目前一些工作［4，18］也考慮到了不同個(gè)體的樣本背后的因果關(guān)系可能存在一定共性，并嘗試從這種混雜樣本中將樣本劃分成不同類別并學(xué)習(xí)因果關(guān)系。文獻(xiàn)［4］與本文工作同樣是面向多元時(shí)間序列數(shù)據(jù)，但其需要指定群體個(gè)數(shù)，并在此基礎(chǔ)上學(xué)習(xí)個(gè)體特定的因果關(guān)系和共性的因果關(guān)系，再基于共性的因果關(guān)系計(jì)算某個(gè)個(gè)體屬于不同群體的概率。文獻(xiàn)［4］與本文工作的不同之處在于，其學(xué)習(xí)的是個(gè)體個(gè)性的因果關(guān)系以及個(gè)體間共性的因果關(guān)系，認(rèn)為不同個(gè)體背后的因果關(guān)系仍然是不同的，僅是存在一定的共性因果關(guān)系，并基于共性進(jìn)行聚類。本文工作則是在個(gè)體間的因果關(guān)系可能相同的場(chǎng)景下，對(duì)多個(gè)個(gè)體聚類并學(xué)習(xí)群體因果關(guān)系。此外，在實(shí)際應(yīng)用中，群體個(gè)數(shù)往往是未知的，而本文所提出的算法無(wú)須指定群體個(gè)數(shù)。文獻(xiàn)［18］考慮的則是二元變量之間的因果關(guān)系，在為每個(gè)樣本識(shí)別了因果關(guān)系后，基于每個(gè)樣本所對(duì)應(yīng)的因果關(guān)系參數(shù)進(jìn)行K-Means［19］聚類，因此也面臨著需要指定群體個(gè)數(shù)的問題。

2 問題定義

本節(jié)對(duì)所研究的問題進(jìn)行符號(hào)化定義和說明。定義n個(gè)個(gè)體的多元時(shí)間序列數(shù)據(jù)集X={X1，X2，…，Xn}，其中第s個(gè)個(gè)體的多元時(shí)間序列數(shù)據(jù)Xs={，，…，}，且每個(gè)個(gè)體的變量 數(shù)均為m，時(shí)間序列長(zhǎng)度均為T。將最長(zhǎng)的因果關(guān)系時(shí)間間隔記為k，多元時(shí)間序列數(shù)據(jù)的因果關(guān)系表示為k+1個(gè)m×m的矩陣{B0，B1，…，Bk}，Bτ[i，j]≠0 表示第t時(shí)刻的變量xi(t)受到第t-τ時(shí)刻的變量xj(t-τ)的因果影響，且時(shí)間間隔為τ，τ∈{0，1，…，k}，其中τ=0時(shí)Bτ表示瞬時(shí)因果關(guān)系，τ＞0時(shí)Bτ表示時(shí)延因果關(guān)系。

本文所考慮的問題是：給定多元時(shí)間序列數(shù)據(jù)集X，如何基于數(shù)據(jù)背后的因果關(guān)系將n個(gè)個(gè)體劃分成c個(gè)群體，且c無(wú)須人工先驗(yàn)指定，并學(xué)習(xí)每個(gè)群體的因果關(guān)系。

3 LEAD 算法

針對(duì)所研究問題，本文提出面向多元時(shí)間序列的群體因果關(guān)系發(fā)現(xiàn)算法LEAD。該算法分為2 個(gè)階段：第一階段是基于因果關(guān)系對(duì)全體的多元時(shí)間序列數(shù)據(jù)進(jìn)行個(gè)體級(jí)別的聚類，將具有相同因果關(guān)系的個(gè)體聚集成一個(gè)群體，并且無(wú)須指定群體的個(gè)數(shù)；第二階段基于前一階段的聚類結(jié)果，利用每個(gè)群體中所有個(gè)體的多元時(shí)間序列數(shù)據(jù)進(jìn)行群體因果關(guān)系發(fā)現(xiàn)。圖1 展示了LEAD 算法的整體框架。

圖1 LEAD 算法框架Fig.1 Framework of LEAD algorithm

下面將先討論如何基于變分推斷的方法，利用個(gè)體的多元時(shí)間序列數(shù)據(jù)學(xué)習(xí)該個(gè)體背后的因果關(guān)系，該方法在LEAD 算法的2 個(gè)階段中都會(huì)涉及，隨后將對(duì)2 個(gè)階段的具體流程進(jìn)行介紹。

3.1 基于變分推斷的個(gè)體因果關(guān)系發(fā)現(xiàn)算法

本文主要考慮變量之間的因果關(guān)系服從線性非高斯的情況，并采用函數(shù)模型進(jìn)行建模。第s個(gè)個(gè)體的多元時(shí)間序列數(shù)據(jù)的生成機(jī)制如下：

其中：N(·)表示一個(gè)高斯分布。因此，主要的推斷任務(wù)是估計(jì)真實(shí)后驗(yàn)分布p(Bs|Xs)。然而直接計(jì)算該后驗(yàn)分布十分困難，本小節(jié)采用變分推斷的技巧，利用一個(gè)易于計(jì)算的變分分布q(Bs)去近似真實(shí)后驗(yàn)分布p(Bs|Xs)。首先給出第s個(gè)個(gè)體的觀測(cè)數(shù)據(jù)Xs的對(duì)數(shù)似然度：

從上述推導(dǎo)中可以看出，觀測(cè)數(shù)據(jù)Xs的對(duì)數(shù)似然度存在一個(gè)下界，當(dāng)且僅當(dāng)KL(q(Bs)||p(Bs|Xs))=0時(shí)，該下界等式成立。因此，只需要最大化這個(gè)下界，就可以找到一個(gè)變分分布q(Bs)足夠近似真實(shí)后驗(yàn)分布p(Bs|Xs)，由此得到以下目標(biāo)函數(shù)：

通過優(yōu)化變分分布q(Bs)來(lái)最大化式（4）即可近似真實(shí)后驗(yàn)分布，從而得到第s個(gè)個(gè)體背后因果關(guān)系。

為了計(jì)算式（4），對(duì)其進(jìn)行分解如下：

將等式最后的第一項(xiàng)記為L(zhǎng)ell(q(Bs)，Xs)，第二項(xiàng)記為L(zhǎng)kl(q，p)，即：

可以看出：Lell(q(Bs)，Xs)是給定變分分布q(Bs)下觀測(cè)數(shù)據(jù)Xs的對(duì)數(shù)似然度的期望；Lkl(q，p)是近似后驗(yàn)分布q(Bs)與先驗(yàn)分布p(Bs)的KL 散度的負(fù)數(shù)。因此，最大化式（4）意味著最大化觀測(cè)數(shù)據(jù)Xs的對(duì)數(shù)似然度且同時(shí)最小化q(Bs)與先驗(yàn)分布的差距。

式（4）的計(jì)算可以通過式（6）和式（7）實(shí)現(xiàn)，因此，對(duì)變分分布q(Bs)做變分推斷中常見的平均場(chǎng)假設(shè)，即：

下面給出lnp(Xs|Bs)的具體計(jì)算方式。由數(shù)據(jù)的產(chǎn)生機(jī)制（式（1））可知：

進(jìn)一步用混合高斯分布來(lái)對(duì)噪聲項(xiàng)進(jìn)行建模，即：

對(duì)于式（7），先驗(yàn)分布p(Bs)和變分分布q(Bs)都是高斯分布，因此，可以通過解析式直接計(jì)算2 個(gè)高斯分布的KL 散度?；谏鲜鐾茖?dǎo)分析，可以對(duì)目標(biāo)函數(shù)式（4）進(jìn)行估計(jì)并計(jì)算梯度，進(jìn)而采用梯度上升的方式對(duì)變分分布q(Bs)和噪聲分布pe(es(t))中的參數(shù)進(jìn)行優(yōu)化，從而最大化目標(biāo)函數(shù)，得到第s個(gè)個(gè)體背后的因果關(guān)系Bs。

3.2 基于因果關(guān)系的聚類算法

上節(jié)給出了如何基于變分推斷的方法，利用個(gè)體的多元時(shí)間序列數(shù)據(jù)學(xué)習(xí)該個(gè)體背后的因果關(guān)系，本節(jié)將介紹如何基于學(xué)習(xí)到的個(gè)體因果關(guān)系進(jìn)行聚類。

設(shè)計(jì)聚類算法的目的是為了將具有相同或者相似因果關(guān)系的個(gè)體聚成一個(gè)群體，從而利用群體的數(shù)據(jù)來(lái)學(xué)習(xí)群體背后的因果關(guān)系。因此，本文從因果關(guān)系的角度來(lái)度量個(gè)體之間的相似性并設(shè)計(jì)個(gè)體聚類過程，提出基于因果關(guān)系的時(shí)間序列聚類算法，如算法1 所示。

算法1基于因果關(guān)系的時(shí)間序列聚類算法

算法2Split(Cv，Lv)函數(shù)算法

算法1 的設(shè)計(jì)主要借鑒了分裂層次聚類的思想，初始時(shí)將全部個(gè)體作為一個(gè)群體，隨機(jī)選擇一個(gè)個(gè)體Xs作為群體代表進(jìn)行因果關(guān)系學(xué)習(xí)，以學(xué)習(xí)得到的變分分布q(Bs)作為群體共同的因果關(guān)系，并計(jì)算每個(gè)個(gè)體在群體因果關(guān)系下的似然度（算法第1～7 行）。個(gè)體似然度高意味著該個(gè)體的因果關(guān)系與該個(gè)體所在群體對(duì)應(yīng)的群體共同因果關(guān)系相似；個(gè)體似然度低意味著該個(gè)體的因果關(guān)系與該個(gè)體所在群體對(duì)應(yīng)的群體共同因果關(guān)系不相似。隨后算法進(jìn)入循環(huán)過程（算法第8～18 行），在每一次循環(huán)中，需要從當(dāng)前已有的群體中選擇出群體內(nèi)部個(gè)體似然度方差最大的群體。群體內(nèi)部個(gè)體似然度的方差越大，意味著該群體內(nèi)部個(gè)體的數(shù)據(jù)背后越可能是不同的因果關(guān)系，導(dǎo)致在相同的群體因果關(guān)系下，有些個(gè)體似然度高，有些個(gè)體似然度低。在此基礎(chǔ)上，基于選出的群體，利用算法2 對(duì)其做分裂，并根據(jù)劃分前后群體所有個(gè)體的似然度之和是否增加來(lái)決定是否接受這次分裂，如果總體似然度之和沒有增加，則拒絕分裂且跳出循環(huán)，此時(shí)算法1 返回聚類結(jié)果并結(jié)束（算法第19 和20 行）。

由于分裂后的2 個(gè)群體分別學(xué)習(xí)各自的因果關(guān)系并在對(duì)應(yīng)群體因果關(guān)系下計(jì)算該群體內(nèi)的個(gè)體的似然度，無(wú)法直接比較分裂前后群體方差的變化，同時(shí)本文希望個(gè)體似然度高，因此算法1 以群體內(nèi)部的所有個(gè)體似然度之和是否增加作為是否分裂的條件。同時(shí)，算法1 在拒絕分裂后隨即跳出循環(huán)結(jié)束算法流程具有2 點(diǎn)好處：一方面緩解了經(jīng)典分裂層次聚類算法需要分裂到單個(gè)個(gè)體作為一個(gè)簇才停止的情況，降低了時(shí)間復(fù)雜度；另一方面無(wú)須為算法指定群體的個(gè)數(shù)，循環(huán)停止時(shí)即得到了群體的個(gè)數(shù)，避免了經(jīng)典層次聚類算法結(jié)束后仍需要選擇簇個(gè)數(shù)的問題。

3.3 群體因果關(guān)系發(fā)現(xiàn)算法

上節(jié)基于因果關(guān)系將全部個(gè)體劃分成了多個(gè)群體，并且無(wú)須指定群體的個(gè)數(shù)，本小節(jié)將基于前一階段的聚類結(jié)果，利用每個(gè)群體中的所有個(gè)體的多元時(shí)間序列數(shù)據(jù)進(jìn)行群體因果關(guān)系發(fā)現(xiàn)。

本節(jié)所提出的群體因果發(fā)現(xiàn)算法主要基于貝葉斯思想，即利用數(shù)據(jù)更新先驗(yàn)分布得到后驗(yàn)分布并循環(huán)這一更新過程，算法偽代碼如算法3 所示。

算法3群體因果關(guān)系發(fā)現(xiàn)算法

算法3對(duì)于每個(gè)群體進(jìn)行遍歷，在同一個(gè)群體內(nèi)，先初始化該群體因果關(guān)系的先驗(yàn)分布（算法第2 和3 行），再對(duì)該群體內(nèi)的個(gè)體進(jìn)行遍歷，對(duì)于一個(gè)個(gè)體，利用該個(gè)體數(shù)據(jù)對(duì)式（4）最大化得到近似后驗(yàn)分布，并以此近似后驗(yàn)分布作為下一次的先驗(yàn)分布，重復(fù)這個(gè)過程直到所有個(gè)體都參與過群體因果關(guān)系分布的更新（算法第4 和7 行），以最終得到的近似后驗(yàn)分布的期望作為該群體的群體因果關(guān)系并保存結(jié)果（算法第8～9 行），直到遍歷完所有群體。通過上述步驟就可以學(xué)習(xí)出各個(gè)群體共享的群體因果關(guān)系。得益于3.1 節(jié)中基于變分推斷的個(gè)體因果關(guān)系發(fā)現(xiàn)算法與貝葉斯思想，算法3 能夠充分地利用群體內(nèi)所有個(gè)體的數(shù)據(jù)，并且學(xué)習(xí)出群體背后的因果關(guān)系。

4 實(shí)驗(yàn)

為了驗(yàn)證本文提出的LEAD 算法，本節(jié)將使用仿真數(shù)據(jù)和真實(shí)數(shù)據(jù)對(duì)算法進(jìn)行實(shí)驗(yàn)評(píng)估。

4.1 仿真數(shù)據(jù)集

設(shè)計(jì)4 組控制變量實(shí)驗(yàn)，具體如下：

1）群體個(gè)數(shù)c={2，3，4，5，6}，個(gè)體總數(shù)n=60，個(gè)體的時(shí)間序列長(zhǎng)度T=60，變量個(gè)數(shù)m=5。

2）群體個(gè)數(shù)c=2，個(gè)體總數(shù)n={20，30，40，50，60}，個(gè)體的時(shí)間序列長(zhǎng)度T=60，變量個(gè)數(shù)m=5。

3）群體個(gè)數(shù)c=2，個(gè)體總數(shù)n=30，個(gè)體的時(shí)間序列長(zhǎng)度T={40，60，80，100，120}，變量個(gè)數(shù)m=5。

4）群體個(gè)數(shù)c=2，個(gè)體總數(shù)n=30，個(gè)體的時(shí)間序列長(zhǎng)度T=60，變量個(gè)數(shù)m={6，8，10，12，14}。

本文所提出的LEAD 算法對(duì)多個(gè)個(gè)體進(jìn)行聚類并學(xué)習(xí)個(gè)體背后的因果關(guān)系（即該個(gè)體所在群體的群體因果關(guān)系），因此，本節(jié)將其與主流聚類算法對(duì)比聚類的效果，與時(shí)間序列因果發(fā)現(xiàn)算法對(duì)比因果發(fā)現(xiàn)的效果。

聚類對(duì)比算法如下：1）K-Means（Euclidean），基于歐氏距離的K-Means 算法［19］；2）K-Means（DTW），基于DTW 距離［22］的K-Means 算法；3）DBSCAN（DTW），基于DTW 距離的DBSCAN 算法［23］；4）OPTICS（DTW），基于DTW 距離的OPTICS 算法［24］。上述聚類對(duì)比算法采用Python 包scikit-learn 中的實(shí)現(xiàn)方法，參數(shù)設(shè)置為各算法實(shí)現(xiàn)的默認(rèn)參數(shù)。其中，K-Means 算法需要指定簇的個(gè)數(shù)K，本小節(jié)實(shí)驗(yàn)將按真實(shí)值來(lái)指定K，注意K在實(shí)際中往往是未知的。

因果關(guān)系發(fā)現(xiàn)對(duì)比算法如下：1）基于約束的PCMCI 算法［8］；2）基于評(píng)分的DYNOTEARS 算法［9］；3）基于函數(shù)的VAR-LiNGAM 算法［10］。上述時(shí)間序列因果發(fā)現(xiàn)對(duì)比算法均采用官方開源代碼實(shí)現(xiàn)，并按真實(shí)值設(shè)置對(duì)比算法的最大因果關(guān)系時(shí)間間隔參數(shù)，其余參數(shù)采用原算法實(shí)現(xiàn)的默認(rèn)參數(shù)。

在評(píng)價(jià)指標(biāo)方面，本實(shí)驗(yàn)選擇ARI（Adjusted Rand Index）指標(biāo)［25］對(duì)聚類算法進(jìn)行評(píng)價(jià)，選擇AUC（Area Under Curve）指標(biāo)對(duì)因果關(guān)系發(fā)現(xiàn)算法進(jìn)行評(píng)價(jià)。ARI 計(jì)算公式為：

其中：n是給定數(shù)據(jù)集中個(gè)體總數(shù)；nij是在聚類結(jié)果R中屬于簇i同時(shí)在聚類結(jié)果R*中屬于簇j的個(gè)體數(shù)量；ni是在聚類結(jié)果R中屬于簇i的個(gè)體數(shù)量；nj是在聚類結(jié)果R*中屬于簇j的個(gè)體數(shù)量；

ARI 的取值范圍是［-1，1］，隨機(jī)劃分下ARI 為0，ARI 越大說明R和R*結(jié)果越一致。AUC 是ROC曲線下的面積，ROC 曲線是以假正率（FPR）為橫軸，真正率（TPR）為縱軸，在多個(gè)閾值下得到一系列點(diǎn)（FPR，TPR）所構(gòu)成的曲線，其中FPR 和TPR 計(jì)算公式為：

其中：將預(yù)測(cè)因果結(jié)構(gòu)與實(shí)際因果結(jié)構(gòu)相比；FP是預(yù)測(cè)有邊但實(shí)際無(wú)邊的數(shù)量；TN是預(yù)測(cè)無(wú)邊實(shí)際也無(wú)邊的數(shù)量；TP是預(yù)測(cè)有邊實(shí)際也有邊的數(shù)量；FN是預(yù)測(cè)無(wú)邊但實(shí)際有邊的數(shù)量，AUC 的取值范圍是［0，1］，AUC 越接近1，預(yù)測(cè)因果結(jié)構(gòu)與真實(shí)因果結(jié)構(gòu)越接近。

在仿真實(shí)驗(yàn)中，相同參數(shù)的實(shí)驗(yàn)都在10 個(gè)不同隨機(jī)種子下進(jìn)行數(shù)據(jù)生成并運(yùn)行LEAD 算法和對(duì)比算法，分別計(jì)算評(píng)價(jià)指標(biāo)并取10 次結(jié)果的平均值作為算法的最終結(jié)果，結(jié)果如圖2～圖5 所示。

圖2 針對(duì)群體個(gè)數(shù)的控制實(shí)驗(yàn)Fig.2 Control experiments on the number of groups

圖2 顯示了個(gè)體總數(shù)相同但背后群體個(gè)數(shù)不同情況下的實(shí)驗(yàn)結(jié)果?？梢钥闯霰疚奶岢龅腖EAD 算法的因果發(fā)現(xiàn)效果要優(yōu)于對(duì)比算法。但隨著群體個(gè)數(shù)的增加，LEAD 算法的因果發(fā)現(xiàn)效果下降，這是因?yàn)閷?shí)際群體個(gè)數(shù)越多，聚類的難度會(huì)越大，特別是在群體個(gè)數(shù)未知的情況下難度會(huì)更大，因此LEAD 算法第一階段的聚類效果變差，導(dǎo)致第二階段的因果發(fā)現(xiàn)效果也下降。事實(shí)上，圖2（b）中的聚類方法都隨著群體個(gè)數(shù)的增加而效果變差，其中K-Means（Euclidean）和K-Means（DTW）由于給定了真實(shí)群體個(gè)數(shù)作為算法參數(shù)，因此效果要優(yōu)于其他方法。而本文提出的LEAD 算法在沒有給定真實(shí)群體個(gè)數(shù)的情況下，也取得了相近的聚類效果。

圖3 顯示了群體個(gè)數(shù)相同但個(gè)體總數(shù)不同情況下的實(shí)驗(yàn)結(jié)果?？梢钥闯鯨EAD 算法的因果發(fā)現(xiàn)能力遠(yuǎn)高于對(duì)比算法，這是因?yàn)長(zhǎng)EAD 算法可以利用群體內(nèi)多個(gè)個(gè)體數(shù)據(jù)。隨著個(gè)體總數(shù)增加，聚類的難度增大，LEAD 算法的聚類效果有所影響，但其因果發(fā)現(xiàn)效果基本沒有影響，甚至因?yàn)閿?shù)據(jù)增加而略微改善。

圖3 針對(duì)個(gè)體總數(shù)的控制實(shí)驗(yàn)Fig.3 Control experiments on the total number of individuals

圖4 顯示了不同時(shí)間序列長(zhǎng)度下的實(shí)驗(yàn)結(jié)果。隨著序列長(zhǎng)度增加，圖4 中的因果發(fā)現(xiàn)算法和聚類算法效果都有所改善，且LEAD 算法無(wú)論是因果發(fā)現(xiàn)還是聚類的效果均要優(yōu)于對(duì)比算法，說明LEAD算法不僅具有優(yōu)異的因果發(fā)現(xiàn)能力，而且還具有很好的多元時(shí)間序列聚類能力。

圖4 針對(duì)時(shí)間序列長(zhǎng)度的控制實(shí)驗(yàn)Fig.4 Control experiments on the length of time series

圖5顯示了不同變量個(gè)數(shù)下的實(shí)驗(yàn)結(jié)果。從圖5（b）可以看出，隨著變量個(gè)數(shù)增加，個(gè)體間的差異信息越來(lái)越豐富，聚類的難度會(huì)下降，聚類算法的效果均明顯改善。但與此同時(shí)，變量間的因果關(guān)系也更復(fù)雜，因果發(fā)現(xiàn)的難度增大，因此圖5（a）中的因果發(fā)現(xiàn)算法均隨著變量個(gè)數(shù)增加而效果變差，但LEAD 算法仍然優(yōu)于對(duì)比算法，驗(yàn)證了算法的有效性。

圖5 針對(duì)變量個(gè)數(shù)的控制實(shí)驗(yàn)Fig.5 Control experiments on the number of variables

4.2 真實(shí)數(shù)據(jù)集

本節(jié)選擇在真實(shí)的Sachs 數(shù)據(jù)集［5］進(jìn)行測(cè)試。Sachs 數(shù)據(jù)集在不同的干預(yù)措施下測(cè)量了細(xì)胞中11種磷酸化蛋白質(zhì)和磷脂分子的濃度變化，不同干預(yù)措施下這11 種分子的因果關(guān)系會(huì)發(fā)生變化［4］，本節(jié)使用干預(yù)cd3cd28+U0126 和干預(yù)cd3cd28+aktinhib下的數(shù)據(jù)，2 種干預(yù)下的數(shù)據(jù)分別屬于2 個(gè)群體，并以30 作為一個(gè)個(gè)體的時(shí)間序列長(zhǎng)度對(duì)數(shù)據(jù)進(jìn)行切分，得到每個(gè)25 個(gè)個(gè)體的數(shù)據(jù)。

將LEAD 算法和聚類算法應(yīng)用到該數(shù)據(jù)集中，并展示LEAD 算法所學(xué)習(xí)出的因果關(guān)系。由表1 可以看出，本文提出的LEAD 算法的聚類效果要優(yōu)于對(duì)比算法，因?yàn)長(zhǎng)EAD 算法能很好地捕捉到因果關(guān)系上的變化而其他距離度量方式不容易捕捉到這種變化。進(jìn)一步分析聚類結(jié)果發(fā)現(xiàn)，LEAD 算法很好地將屬于干預(yù)cd3cd28+aktinhib 組的25 個(gè)個(gè)體劃分在同一群體（記為群體a1）中，但是將屬于干預(yù)cd3cd28+U0126 組的25 個(gè)個(gè)體分成了2 個(gè)群體，其中一個(gè)群體（記為群體U1）有21 個(gè)個(gè)體，另一個(gè)群體（記為群體U2）有4 個(gè)個(gè)體。

表1 在真實(shí)數(shù)據(jù)集中的聚類性能 Table 1 The clustering performance on real dataset

本文算法學(xué)習(xí)的因果關(guān)系如圖6 所示，其中有向邊a1 為群體a1 中的因果關(guān)系，有向邊U1 為群體U1 中的因果關(guān)系，有向邊U2 為群體U2 中的因果關(guān)系，僅顯示每個(gè)群體中影響權(quán)重前十的因果關(guān)系?？梢钥闯?，群體a1 的重要因果關(guān)系與群體U1、群體U2有較大區(qū)別，而群體U1、U2的重要因果關(guān)系有5條相同，說明算法能夠區(qū)分在不同干預(yù)下的個(gè)體。此外，在群體U1、U2 中都發(fā)現(xiàn)了因果關(guān)系Erk →Akt，該因果關(guān)系與文獻(xiàn)［5］中報(bào)告的一致，但在群體a1 中并未被發(fā)現(xiàn)，這是因?yàn)樵谌后wa1 中Akt 被干預(yù)，從而切斷了Akt 的原因變量，驗(yàn)證了算法的有效性。

圖6 真實(shí)數(shù)據(jù)集中學(xué)習(xí)到的3 個(gè)群體因果關(guān)系Fig.6 Three collective causal relations learned from real dataset

5 結(jié)束語(yǔ)

現(xiàn)有因果關(guān)系發(fā)現(xiàn)方法在多個(gè)個(gè)體數(shù)據(jù)背后有相同因果關(guān)系的情況下樣本利用不足。針對(duì)該問題，本文提出一種面向多元時(shí)間序列的群體因果發(fā)現(xiàn)算法。該算法通過基于因果關(guān)系的聚類，將多個(gè)個(gè)體按照因果關(guān)系的相同或不同劃分成多個(gè)群體，且算法無(wú)須指定群體個(gè)數(shù)。在此基礎(chǔ)上，為每個(gè)群體通過群體因果發(fā)現(xiàn)算法學(xué)習(xí)群體因果關(guān)系。實(shí)驗(yàn)結(jié)果表明，該算法可以充分利用具有相同因果關(guān)系的多個(gè)個(gè)體數(shù)據(jù)，因果發(fā)現(xiàn)能力優(yōu)于對(duì)比算法，并且同時(shí)具有和對(duì)比算法相近的多元時(shí)間序列聚類能力。本文關(guān)注于因果關(guān)系是線性的情況，下一步將考慮在因果關(guān)系滿足非線性的情況下如何進(jìn)行群體因果發(fā)現(xiàn)，進(jìn)一步提升算法的適用范圍。