曾 艷，郝志峰,2+，蔡瑞初，謝 峰

(1.廣東工業(yè)大學(xué) 計算機學(xué)院，廣東 廣州 510006；2.佛山科學(xué)技術(shù)學(xué)院 數(shù)學(xué)與大數(shù)據(jù)學(xué)院， 廣東 佛山 528000；3.北京大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，北京 100080)

0 引 言

傳統(tǒng)的因果關(guān)系發(fā)現(xiàn)算法多數(shù)專注于探索可觀測變量之間的因果關(guān)系[1-6]，然而，在現(xiàn)實生活中，事物本質(zhì)的因果關(guān)系并不是位于可觀測變量中，而是位于隱變量(latent variables)中，例如焦慮、失望或處理事物的能力等[2,7]。因此，從觀察數(shù)據(jù)中發(fā)現(xiàn)隱變量間的因果關(guān)系近年來獲得了不少研究學(xué)者的興趣[7-9]。

經(jīng)典的從觀測變量中發(fā)現(xiàn)隱變量的因果關(guān)系算法可以分為兩類：基于高斯性和基于非高斯性的算法?；诟咚剐缘乃惴╗7,10]主要利用了數(shù)據(jù)結(jié)構(gòu)中的協(xié)方差信息，進而確定隱變量之間的馬爾科夫等價類結(jié)構(gòu)?；诜歉咚沟乃惴╗9,11]等使用了數(shù)據(jù)的非高斯性，解決了該不可唯一識別的問題。這些工作專注于研究靜態(tài)貝葉斯網(wǎng)絡(luò)的模型。然而，在實際應(yīng)用領(lǐng)域中，特別是在經(jīng)濟學(xué)中，因果效應(yīng)往往在動態(tài)貝葉斯網(wǎng)絡(luò)中體現(xiàn)，而不僅僅發(fā)生在瞬時的因果關(guān)系中。如果此時忽略考慮瞬時效應(yīng)或者延時效應(yīng)，那么在因果結(jié)構(gòu)的估計中將會引入額外的錯誤而導(dǎo)致不合理的解釋。現(xiàn)有的基于結(jié)構(gòu)向量自回歸模型(structural vector autoregression model，SVAR)[12]的工作考慮了可觀測變量的瞬時性和延時性的因果效應(yīng)，但沒有考慮隱變量的情況。

因此，本文首先構(gòu)建了時序隱變量模型(latent factor models for time-series data，LFTS)，考慮了隱變量之間的瞬時性和延時性的因果效應(yīng)，提出了基于LFTS模型的因果關(guān)系發(fā)現(xiàn)算法。該算法采用了三階段的方式估計模型的網(wǎng)絡(luò)結(jié)構(gòu)，包括隱變量之間的瞬時因果效應(yīng)矩陣與延時因果效應(yīng)矩陣。最后，通過仿真數(shù)據(jù)實驗，表明了算法的有效性。

1 相關(guān)工作

1.1 結(jié)構(gòu)向量自回歸模型

自回歸模型(autoregressive model，AR)[12]是統(tǒng)計上一種處理時間序列數(shù)據(jù)的方法，被廣泛運用在經(jīng)濟學(xué)、資訊學(xué)、自然現(xiàn)象等領(lǐng)域上。結(jié)構(gòu)向量自回歸模型(structural vector autoregression model，SVAR)是AR模型由一個變量到多個變量的擴展，同時它不僅考慮了各個變量之間的延時效應(yīng)，而且考慮了即時性的結(jié)構(gòu)關(guān)系。值得注意的是，傳統(tǒng)的AR模型和貝葉斯網(wǎng)絡(luò)(Bayesian networks，BN)或結(jié)構(gòu)方程模型(structural equation models，SEM)均屬于SVAR模型的特殊形式，SVAR模型是AR模型與SEM模型的結(jié)合體。SVAR模型旨在學(xué)習(xí)隨機變量之間的延時效應(yīng)矩陣與瞬時效應(yīng)矩陣，從而達到預(yù)測的目的或者實現(xiàn)因果關(guān)系的分析，提高對數(shù)據(jù)的可解釋性。

其數(shù)學(xué)表示為：假定在一個平穩(wěn)時間序列上有M個獨立的樣本，xm,t表示第t個時間序列的第m個數(shù)據(jù)樣本，則其p階SVAR模型可以如式(1)所示

xm,t=B0xm,t+B1xm,t-1+…+Bpxm,t-p+em,t

(1)

其中，t∈{p,…,T}，m∈{1,…,M}，em,t表示對應(yīng)的噪聲變量且它們之間在不同時間序列和在同一個時間序列上互相獨立。Bi,i∈{1,…,p} 表示延時效應(yīng)矩陣，而B0表示瞬時效應(yīng)矩陣，它們屬于SVAR模型的回歸參數(shù)。

1.2 隱變量模型

隱變量模型[7]，也稱為隱因子模型(latent variable/factor models，LVM/LFM)，是統(tǒng)計學(xué)上用以描述可觀測變量和不可觀測變量(隱變量/隱因子)之間關(guān)系的經(jīng)典模型。該模型被廣泛應(yīng)用于心理學(xué)、社會學(xué)和教育研究等領(lǐng)域中。隱變量模型通常假設(shè)可觀測變量由隱因子產(chǎn)生。特別地，當(dāng)所有變量都是連續(xù)的隨機變量時，從SEM的角度上看，隱變量模型可以表示為如下的數(shù)學(xué)表達式

fm=Bfm+em

(2)

xm=Gfm+nm

(3)

其中，fm,xm分別為隱變量向量與可觀測變量向量的第m個樣本；em,nm分別為fm與xm的噪音變量向量的第m個樣本。B是表示連接強度矩陣而G是因子載荷矩陣。B表示隱變量之間的因果結(jié)構(gòu)關(guān)系，而G表示可觀測變量與隱變量之間的聯(lián)系。通常地，式(2)被稱為該隱變量模型中的結(jié)構(gòu)模型，而式(3)被稱為測量模型。

由于隱變量的不可觀測性，該模型的估計一直以來被認為是一個難點和挑戰(zhàn)。早前的估計工作主要利用了數(shù)據(jù)的高斯性，然而這類方法會產(chǎn)生許多不可識別的結(jié)構(gòu)模型[7,10]；為了克服不可識別性的問題，近年來，學(xué)者們專注于利用數(shù)據(jù)的非高斯性，實現(xiàn)模型的唯一識別性[8,9]。然而，上述方法主要估計的是靜態(tài)的隱變量模型，而具有時序信息的動態(tài)數(shù)據(jù)在現(xiàn)實生活中卻是常見的。因此，研究時序隱變量模型的因果關(guān)系發(fā)現(xiàn)算法是必要且重要的。

2 基于時序隱因子模型的因果關(guān)系發(fā)現(xiàn)算法

本節(jié)在相關(guān)工作的基礎(chǔ)上建立了面向時序數(shù)據(jù)的隱因子模型，接著提出了基于該模型的因果關(guān)系發(fā)現(xiàn)算法。

2.1 時序隱變量模型

本文將隱變量模型擴展到時間序列上，提出了時序隱變量模型，簡稱為LFTS，其數(shù)學(xué)表達式具體如下

fm,t=B0fm,t+B1fm,t-1+…+Bpfm,t-p+em,t

(4)

xm,t=Gfm,t+nm,t

(5)

其中，fm,t,xm,t表示第t個時間序列隱變量、可觀測變量的第m個數(shù)據(jù)樣本，t∈{p,…,T}，m∈{1,…,M}， 自回歸階數(shù)為p階，Bi,i∈{1,…,p} 表示隱變量的延時效應(yīng)矩陣，而B0表示隱變量之間的瞬時效應(yīng)矩陣。值得注意的是，在t時刻，隱變量之間的因果關(guān)系網(wǎng)絡(luò)是一個有向無環(huán)圖，所以在經(jīng)過一系列的矩陣變換操作之后，瞬時效應(yīng)矩陣B0可以變換成為一個嚴(yán)格的下三角矩陣。em,t,nm,t分別為隱變量、可觀測變量的噪音變量向量在第t個時間序列的第m個樣本，G是因子載荷矩陣。類似地，本文稱式(4)為LFTS的結(jié)構(gòu)模型，而式(5)為其測量模型。

為了獲得模型的識別性[9,13]，假設(shè)噪音向量em,t,nm,t在時刻t內(nèi)或者在時間之間是互相獨立的；噪音變量服從非高斯分布；每一個隱變量至少產(chǎn)生兩個純測量變量(純測量變量指的是只由一個隱變量父親產(chǎn)生的可觀測變量[2])。式(4)可以重寫為如下向量自回歸模型的形式

(I-B0)fm,t=B1fm,t-1+…+I-Bpfm,t-p+em,t

(6)

fm,t=A1fm,t-1+…+Apfm,t-p+εm,t

(7)

其中

Ai=(I-B0)-1Bi

(8)

εm,t=(I-B0)-1em,t

(9)

I為單位矩陣。由此可知，式(7)是向量自回歸模型的簡化型。除此之外，由式(9)可得

εm,t=B0εm,t+em,t

(10)

其中

εm,t=fm,t-A1fm,t-1-…-Apfm,t-p

(11)

由式(10)以及本文所做的假設(shè)可知，LFTS模型中的隱變量的殘差εm,t在本質(zhì)上服從線性非高斯無環(huán)模型(linear non-Gaussian acyclic models，LiNGAM)[14]。而式(5)的測量模型直接是隱因子模型的簡化型。

圖1提供了本模型的一個例子，即，一個具有3個隱變量、每個隱變量下具有2個純測量變量、自回歸階數(shù)為2的時序隱變量模型。如圖1所示，不同的底色顯示不同的時間序列，即時間為(t),(t-1) 和 (t-2) 以及p=2。 方形的節(jié)點表示隱變量，圓形的節(jié)點表示可觀測變量。fi或xi表示第i個隱變量或可觀測變量。隱變量之間的實線箭頭表示瞬時效應(yīng)矩陣B0， 而虛線箭頭表示延時效應(yīng)矩陣Bi。 隱變量與可觀測變量之間的關(guān)系使用實線箭頭的G表示。由圖1可知，以第m個樣本為例， (t)時刻下的f1受到 (t-1) 時刻下的f1和f3， 以及(t-2) 時刻下的f1共同影響； (t) 時刻下的f2受到同時刻的f1、 (t-1) 時刻下的f2以及 (t-2) 時刻下的f2共同影響。所以，可以定義問題如下：

圖1 本文時序隱變量模型

2.2 LFTS算法

(12)

式(12)由式(5)推導(dǎo)而成，且在噪聲變量足夠小時或者測量變量的個數(shù)足夠多時成立。

由上一小節(jié)分析可知，式(4)表示的是結(jié)構(gòu)向量自回歸模型，它不僅刻畫了隱變量之間的瞬時效應(yīng)矩陣B0， 而且描述了不同時序下隱變量之間的延時效應(yīng)矩陣Bi。 式(4)可以轉(zhuǎn)換成式(7)。由式(7)可知，隱變量fm,t可以表示為過去p個時間序列隱變量的線性組合形式和噪聲εm,t之和，也就是，在本質(zhì)上式(7)屬于經(jīng)典的向量自回歸模型的簡化型，其中矩陣Ai是自回歸矩陣，p是隱變量延時的時序個數(shù)，即自回歸階數(shù)，而噪聲εm,t是革新過程。值得注意的是，當(dāng)瞬時延時矩陣B0為零矩陣時，式(4)的模型則不需要經(jīng)過變換，其自身已經(jīng)是一個向量自回歸模型。

εm,t=Wem,t

(13)

其中

W=(I-B)-1

(14)

(15)

2.2.4 算法偽代碼

根據(jù)前文對算法流程的描述，現(xiàn)將LFTS算法的偽代碼列出如算法1所示：

算法1：LFTS算法

輸入：具有時序信息的數(shù)據(jù)集X， 自回歸階數(shù)p。

(1)歸一化處理數(shù)據(jù)集X；

3 實驗結(jié)果與分析

為了驗證和分析本文提出算法的有效性，本文在模擬數(shù)據(jù)集上做了驗證實驗與對比實驗。具體來說，本文從兩個方面驗證算法的有效性：3.1考慮隱變量之間因果延時效應(yīng)的驗證實驗，使用熱力圖描述估計的瞬時效應(yīng)矩陣B0、 延時效應(yīng)矩陣Bi， 體現(xiàn)本算法的有效性；3.2不考慮隱變量之間因果延時效應(yīng)的對比實驗，此時，選取經(jīng)典的LiNGAM[14]與Notears[5]方法作為隱變量之間瞬時因果結(jié)構(gòu)的估計方法，這兩類方法沒有考慮隱變量之間的因果延時效應(yīng)。

模擬數(shù)據(jù)的生成主要通過兩個階段來實現(xiàn)：首先，令自回歸系數(shù)為1，先隨機產(chǎn)生瞬時效應(yīng)矩陣B0、 延時效應(yīng)矩陣Bi和隱變量的因果次序(因果關(guān)系的強度在[-0.7,-0.2]∪[0.2,0.7]之間取得)。給定樣本量，噪聲變量em,t的分布為非高斯分布，根據(jù)式(4)隨機產(chǎn)生LFTS模型中的結(jié)構(gòu)模型的隱變量模擬數(shù)據(jù)fm,t。 緊接著，假定每個隱變量下具有2個純的測量變量，隨機產(chǎn)生[-0.7,-0.2]∪[0.2,0.7]之間的因子載荷矩陣G， 噪聲變量nm,t的分布為高斯分布，再根據(jù)式(5)隨機產(chǎn)生LFTS模型中的測量模型的可觀測變量數(shù)據(jù)xm,t。

3.1 考慮延時效應(yīng)的實驗分析

本小節(jié)主要進行了不同樣本量情況下的模擬數(shù)據(jù)實驗：設(shè)置隱變量的個數(shù)為5，可觀測變量個數(shù)為10，樣本量為100,200,500,1000。由于本算法旨在發(fā)現(xiàn)時序隱變量模型的網(wǎng)絡(luò)因果結(jié)構(gòu)關(guān)系及其相應(yīng)的因果強度，本文將真實的B0,Bi和估計的B0,Bi作為輸出結(jié)果作對比，即使用熱力圖描述實驗輸出的瞬時效應(yīng)矩陣B0和延時效應(yīng)矩陣Bi。 圖2～圖5展示了算法在不同樣本量下的真實矩陣B0,Bi與估計矩陣B0,Bi的對比圖，其中，圖(a)～圖(d)分別為真實的B0、 估計的B0、 真實的B1和估計的B1的熱力圖。當(dāng)真實的因果效應(yīng)矩陣圖(a)、圖(c)與估計的效應(yīng)矩陣圖(b)、圖(d)越相似，則實驗效果越好。當(dāng)樣本量為100,200和500時，實驗設(shè)置了不同的隨機種子，所以每次實驗產(chǎn)生的真實結(jié)構(gòu)都會發(fā)生變化。由圖2至圖4所示，總體來看，算法估計的網(wǎng)絡(luò)結(jié)構(gòu)圖(b)和圖(d)的因果強度與真實的因果強度圖(a)和圖(c)是令人滿意的，不管是對于瞬時的隱變量之間的結(jié)構(gòu)B0來說，還是對于延時的不同時序隱變量之間的結(jié)構(gòu)Bi來說。具體地說，當(dāng)因果網(wǎng)絡(luò)結(jié)構(gòu)趨于稠密時，如圖2～圖4中的B1， 算法可能會估計少或多1-2條邊，這可能是因為剪枝算法中的剪枝強度過大或過小。在現(xiàn)實生活中，可以根據(jù)實際情況或數(shù)據(jù)的背景知識，選擇合適的剪枝算法避免估計出缺失或者多余的因果結(jié)構(gòu)。

圖2 樣本量為100的真實矩陣B0,Bi和估計矩陣B0,Bi對比結(jié)果

圖3 樣本量為200的真實矩陣B0,Bi和估計矩陣B0,Bi對比結(jié)果

圖4 樣本量為500的真實矩陣B0,Bi和估計矩陣B0,Bi對比結(jié)果

圖5 樣本量為1000的真實矩陣B0,Bi和估計矩陣B0,Bi對比結(jié)果

此外，當(dāng)樣本量為1000時，設(shè)置了和樣本量為500的實驗一樣的隨機種子，可以發(fā)現(xiàn)，當(dāng)固定隨機種子時，隨著樣本量的增大，實驗的效果有一定的提升。特別是對于延時效應(yīng)矩陣Bi來說，其學(xué)習(xí)到的因果結(jié)構(gòu)更加準(zhǔn)確，且因果強度更貼近于真實值。所以，總體來說，算法在時序隱變量模型的因果關(guān)系發(fā)現(xiàn)問題上具有較好的實驗效果。

3.2 不考慮延時效應(yīng)的實驗分析

為了說明考慮因果延時效應(yīng)的必要性以及本算法在考慮因果延時效應(yīng)方面的有效性，本小節(jié)選取了經(jīng)典的LiNGAM[14]與Notears[5]方法作為隱變量之間瞬時因果結(jié)構(gòu)的估計方法，這兩類方法均沒有考慮隱變量之間的因果延時效應(yīng)。本文采用召回率(Recalls)、準(zhǔn)確率(Precisions)和F1得分(F1)作為因果瞬時效應(yīng)矩陣的評估指標(biāo)。F1得分反映了隱變量因果網(wǎng)絡(luò)結(jié)構(gòu)的總體估計優(yōu)劣，具體為

(16)

其中，TP表示正確估計的邊的個數(shù)；FP表示未估計到的邊的個數(shù)；FN表示錯誤估計的邊的個數(shù)。本小節(jié)設(shè)置了在不同樣本量100,200,500,1000下，不同隱變量個數(shù)3,4,5,6的實驗效果。結(jié)果取50次獨立隨機實驗的平均值。結(jié)果如圖6～圖9所示。

圖6 算法在樣本量為100時的實驗效果

圖8 算法在樣本量為500時的實驗效果

圖9 算法在樣本量為1000時的實驗效果

從圖6～圖9的結(jié)果圖分析發(fā)現(xiàn)，隨著樣本量的逐漸增加，所有算法的3個指標(biāo)都呈現(xiàn)不同程度上升的趨勢，且本文提出的LFTS算法均優(yōu)于其余兩個對比算法，說明本算法的有效性。具體地來說，當(dāng)樣本量小于500時，LFTS算法隨著樣本量的增加而效果提高顯著，當(dāng)樣本量為500或1000時，其效果提升的幅度相對穩(wěn)定，這意味著算法在樣本量大于500時對樣本量的魯棒性提高。此外，LFTS算法幾乎在每個不同設(shè)置下的實驗都優(yōu)于其余算法，充分地說明了考慮延時因果效應(yīng)的重要性。如果忽略了該效應(yīng)，那么在估計瞬時因果效應(yīng)矩陣時，會大大地降低估計的準(zhǔn)確性，影響了最終隱變量因果網(wǎng)絡(luò)結(jié)構(gòu)，包括方向的估計，從而影響了實驗數(shù)據(jù)最終的解釋性。

4 結(jié)束語

從觀測變量中發(fā)現(xiàn)隱變量之間因果關(guān)系的算法近年來獲得越來越多學(xué)者的關(guān)注。本文在現(xiàn)有隱變量模型的基礎(chǔ)上，考慮了隱變量之間的瞬時因果效應(yīng)與延時因果效應(yīng)，創(chuàng)建了時序隱變量模型(latent factor models for time-series data，LFTS)，并提出了一種三階段的因果關(guān)系發(fā)現(xiàn)算法，估計了LFTS模型的因果網(wǎng)絡(luò)結(jié)構(gòu)，包括瞬時因果效應(yīng)矩陣與延時因果效應(yīng)矩陣。實驗結(jié)果表明，本文提出的LFTS算法具有較好的學(xué)習(xí)效果，具有較高的估計瞬時因果效應(yīng)與延時因果效應(yīng)的準(zhǔn)確率。如何使用統(tǒng)一的框架實現(xiàn)同時估計瞬時因果效應(yīng)矩陣與延時因果效應(yīng)矩陣，是下一個階段的工作重點。