姓名：劉成龍 ?學(xué)號(hào)：16020199016

轉(zhuǎn)載自：https://www.jiqizhixin.com/articles/2018-11-07-7爆办，有刪節(jié)。

【嵌牛導(dǎo)讀】：使用高斯過程的因果推理：GP CaKe 的基本思路

【嵌牛鼻子】：GP CaKe

【嵌牛提問】：你知道GP CaKe嗎？

【嵌牛正文】：

我們最近開發(fā)出了一種用于時(shí)間序列數(shù)據(jù)中因果推理的全新方法 [Ambrogioni et al., 2017]黎比。我們稱之為「GP CaKe」，即具有因果核的高斯過程（Gaussian Processes?with Causal Kernels）。這種方法不僅縮寫詞（有「蛋糕」的意思）很美味校焦，而且將向量自回歸模型（VAR）那引人注目的特性與動(dòng)態(tài)因果建模（DCM）優(yōu)雅地結(jié)合到了一起。是的统倒，確實(shí)是一舉兩得寨典！

我們最初是為研究心智中的有效連接（即大腦區(qū)域之間的因果交互研究）而開發(fā)了這種方法，但該方法是完全通用型的房匆，也能用在其它任何地方 [1]耸成。我將通過一個(gè)文章系列來解釋 GP CaKe 背后的思想报亩，本文是其中第一篇。在接下來的文章中墓猎，我將通過逐步講解的方式解釋如何使用我們?cè)?GitHub 上提供的代碼捆昏，之后我還將擴(kuò)展介紹這個(gè)模型的延伸方法。

背景：對(duì)多變量時(shí)間序列的分析

這項(xiàng)研究的研究背景是帶有一個(gè)時(shí)間維度的復(fù)雜系統(tǒng)毙沾。在我們的案例中可能是對(duì)多個(gè)大腦區(qū)域的激活狀態(tài)的連續(xù)測(cè)量骗卜，比如通過 EEG、MEG 或 fMRI [2] 等方式左胞，但也可以是連續(xù)的股票交易列表寇仓、天氣現(xiàn)象以及蛋白質(zhì)濃度變化等等。

在統(tǒng)計(jì)學(xué)和機(jī)器學(xué)習(xí)領(lǐng)域烤宙，建模這種復(fù)雜系統(tǒng)的時(shí)間序列的方法主要有兩種：向量自回歸（VAR）[Lütkepohl, 2005] 和動(dòng)態(tài)系統(tǒng)理論（DST）遍烦，后者通常是通過（隨機(jī)）微分方程（SDE）或差分方程 [Abraham & Shaw, 1983] 實(shí)現(xiàn)的。我們將簡(jiǎn)要介紹這兩種方法躺枕，以激發(fā)你的興趣以及說明 GP CaKe 與它們的關(guān)聯(lián)服猪。如果你已經(jīng)很了解這兩種方法了，可以直接跳過這部分拐云，查看后文有關(guān) GP CaKe 的內(nèi)容罢猪。

向量自回歸

向量自回歸（VAR）最基本的形式是預(yù)測(cè)一個(gè)特定變量 x_j(t) 在時(shí)間 t 的值，該變量是另一些變量 x_i 的（隨機(jī)）函數(shù)叉瘩，它們的關(guān)系如下 [3]：

這個(gè)等式的含義如下：

變量 x_j(t) 的信號(hào)取決于該變量從所有其它變量獲得的輸入膳帕。這種依賴關(guān)系的強(qiáng)度由自回歸系數(shù) a_ij(τ) 確定。參數(shù)?τ 是信號(hào) x_j(t) 和 x_i(t) 之間的延后量（lag）薇缅。綜合起來危彩，這表示一個(gè)變量對(duì)另一個(gè)變量的影響可以在 τ=0 時(shí)為零（舉個(gè)例子），然后這種因果影響會(huì)緩慢增大（即 a_ij(τ) 更大）泳桦，只有當(dāng)延后量變得很大時(shí)才會(huì)再次衰減——這意味著在遙遠(yuǎn)過去發(fā)生的某些事情現(xiàn)在已經(jīng)不再重要了汤徽。如果我們繪制出這些系數(shù)與該延后量的函數(shù)關(guān)系，我們就會(huì)得到所謂的脈沖響應(yīng)函數(shù)（impulse response function/IRF）灸撰。最后泻骤，w(t) 描述的是驅(qū)動(dòng)系統(tǒng)的隨機(jī)「創(chuàng)新量（innovations）」或「沖擊量（shocks）」。它們能夠反映 x_j(t) 的內(nèi)部動(dòng)態(tài)梧奢。比如說狱掂，我們國(guó)家的天氣會(huì)受周邊鄰國(guó)天氣狀況的影響（即它們對(duì)我們的氣候有因果影響），但也會(huì)受到我國(guó)內(nèi)部情況的影響亲轨。如果 a_ij(τ)>0趋惨，我們就可以說 x_i 對(duì) x_j 有因果影響（這一思想的實(shí)際實(shí)現(xiàn)將需要某些顯著性測(cè)試）。這意味著因果關(guān)系的時(shí)間概念：一個(gè)變量的過去能為我們提供有關(guān)另一個(gè)變量的未來的信息惦蚊。這種看待因果關(guān)系的角度也被稱為 Wiener-Granger 因果關(guān)系器虾，有時(shí)也被簡(jiǎn)稱為 Granger 因果關(guān)系 [Bressler & Seth, 2011]讯嫂。通過觀察 IRF，我們可以明確地看到 Granger 因果交互的時(shí)間形狀：

圖 1：在三大金融指數(shù)上的 VAR 分析示例兆沙。左圖是每個(gè)變量的時(shí)間序列欧芽，右圖是最大為 10 個(gè)月的延后量的脈沖響應(yīng)函數(shù)（IRF）。注意其中也包含了自響應(yīng)葛圃。

動(dòng)態(tài)系統(tǒng)理論

顧名思義千扔，動(dòng)態(tài)系統(tǒng)理論（DST）也就是建模了系統(tǒng)的動(dòng)態(tài)。以經(jīng)典的 Ornstein-Uhlenbeck 過程為例库正，如下所示：

這描述了一個(gè)隨機(jī)游走過程——隨著時(shí)間推移曲楚，會(huì)逐漸趨近其平均值 μ。

圖 2：五個(gè)全都（漸進(jìn)地）回到同一平均值 μ=0.8 的 Ornstein-Uhlenbeck 過程褥符，但它們各自的噪聲水平和回到均值的速度都各不相同龙誊。

DST 在動(dòng)態(tài)因果建模（DCM）方面也有應(yīng)用 [Friston, 2009]。盡管 DCM 的大多數(shù)實(shí)現(xiàn)都包含一個(gè)專用于 fMRI 的前向模型喷楣，將 DCM 的應(yīng)用限制在了神經(jīng)成像研究領(lǐng)域趟大，但其核心有一個(gè)通用型的微分方程系統(tǒng)：

注意，其中 X=(x_1, ..., x_n)铣焊。此外护昧，A 是一個(gè)包含 x 中變量之間固定交互的矩陣。其作用類似于 VAR 模型的自回歸系數(shù)粗截，但 DCM 中沒有建模延后量。實(shí)際上捣炬，其動(dòng)態(tài)所受的影響是即時(shí)的熊昌。另外的 B 和 C 項(xiàng)表示（節(jié)點(diǎn)特定的）外源輸入 μ，我們這里不會(huì)對(duì)此詳做討論湿酸，因?yàn)?GP CaKe [4] 中（尚）沒有與這些項(xiàng)類同的地方婿屹。

連續(xù)和動(dòng)態(tài)的向量自回歸：GP CaKe

前提內(nèi)容已經(jīng)足夠了，現(xiàn)在該進(jìn)入正題了 [5]推溃。VAR 模型的難點(diǎn)是在實(shí)際操作中昂利，我們沒有足夠多的觀察來可靠地估計(jì)自回歸系數(shù)。由此造成我們的脈沖響應(yīng)函數(shù)充滿噪聲且難以解讀铁坎。此外蜂奸，VAR 模型只能粗略地描述系統(tǒng)的動(dòng)態(tài)。高階交互會(huì)被完全忽視硬萍。DCM 確實(shí)能更廣泛地考慮動(dòng)態(tài)扩所，但是它卻不能建模一個(gè)變量的變化和另一個(gè)變量的動(dòng)態(tài)的變化之間的延遲情況，這是有問題的朴乖。DCM 的某些變體確實(shí)包含一個(gè)延后項(xiàng)祖屏，但會(huì)將其設(shè)置為一個(gè)常量項(xiàng)助赞，而不是我們估計(jì)交互系數(shù)的值的一個(gè)區(qū)間。你可能也已經(jīng)猜到了：GP CaKe 實(shí)際上將延后的交互與動(dòng)態(tài)系統(tǒng)結(jié)合了起來袁勺。下面來一窺究竟雹食。

GP CaKe 的組成可寫成如下形式：

其中，D_j 是微分算子（即其描述了直到第 p 個(gè)導(dǎo)數(shù)的動(dòng)態(tài)）期丰，w_j(t) 仍然是指「創(chuàng)新量（innovation）」或「沖擊量（shock）」項(xiàng)群叶，關(guān)鍵的 C_j(t) 是來自其它變量 i≠j 的因果效應(yīng)的總和。

其中咐汞，C_i(t) 到 C_j(t) 是因果脈沖響應(yīng)函數(shù)（CIRF）盖呼，描述了從 i 到 j 的因果交互。從中可以看到化撕，C_j(t) 是時(shí)間序列的求和（在所有輸入變量上）几晤，并且這些時(shí)間序列與它們的脈沖響應(yīng)函數(shù)進(jìn)行了卷積。這個(gè)定義完全類似于 VAR 模型等式右側(cè)的第一項(xiàng)植阴，但卻是連續(xù)的蟹瘾，而非離散的。然而掠手，GP CaKe 并不簡(jiǎn)單地是 VAR 的連續(xù)式變體憾朴。微分算子 D_j 看似沒啥作用，但實(shí)際非常關(guān)鍵喷鸽。它描述的是一個(gè)變量的內(nèi)部動(dòng)態(tài)众雷，且無論它從其它變量那里得到了什么輸入——而且我們尚未描述這些動(dòng)態(tài)是什么！存在一些（實(shí)際上數(shù)量很有限）可能的選項(xiàng)做祝，比如砾省，這些動(dòng)態(tài)可以是我們之前見過的簡(jiǎn)單的 Ornstein-Uhlenbeck 隨機(jī)游走，即一個(gè)振蕩過程混槐。不管是在什么案例中编兄，都要記住這樣一個(gè)重點(diǎn)：GP CaKe 假設(shè)來自其它變量的輸入會(huì)通過因果脈沖響應(yīng)函數(shù)影響動(dòng)態(tài) D_jx_j(t)，而不只會(huì)影響 x_j(t) 本身声登！

讓我們實(shí)現(xiàn)它

在下一篇文章中狠鸳，我將解釋如何計(jì)算因果脈沖響應(yīng)函數(shù)，這與我們之前的一篇有關(guān)高斯過程回歸中傅立葉變換的文章有關(guān)：https://www.mindcodec.com/the-fourier-transform-through-the-lens-of-gaussian-process-regression/悯嗓。目前件舵，我們就假設(shè)我們已經(jīng)有能幫助我們完成這項(xiàng)任務(wù)的工具了（實(shí)際上我們確實(shí)有，參見對(duì)應(yīng)的 GitHub 項(xiàng)目）脯厨。我們這里只是為了進(jìn)行一點(diǎn)演示芦圾，不會(huì)涉及太多細(xì)節(jié)。我們將使用一個(gè)已知的脈沖響應(yīng)函數(shù)來生成某些數(shù)據(jù)俄认，然后會(huì)嘗試使用 VAR 和 GP CaKe 來恢復(fù)它个少。注意洪乍，目前已經(jīng)存在一些比標(biāo)準(zhǔn)的非正則化的 VAR 更好的實(shí)現(xiàn)，但這里只是用于解釋說明夜焦，普通的 VAR 就足夠了壳澳。

我們從兩個(gè)變量 x_A 和 x_B 開始，它們具有如下的因果關(guān)系：

其中 τ 仍然是兩個(gè)變量之間的時(shí)間延后量茫经，而 s 則是脈沖響應(yīng)的長(zhǎng)度范圍（這個(gè)函數(shù)的形狀如下圖中的紅線所示）巷波。為了我們當(dāng)前的目的，這是一個(gè)任意的變量卸伞，我們只是隨便選取了某個(gè)值抹镊。此外，我們假設(shè)這兩個(gè)變量的內(nèi)部動(dòng)態(tài)都是 Ornstein-Uhlenbeck 過程荤傲，因此

其中垮耳，α 是該過程的弛豫系數(shù)，指示了該時(shí)間序列回到其均值（零）的速度遂黍。我們?yōu)樵搫?dòng)態(tài)系統(tǒng)生成了 100 個(gè)樣本终佛，總長(zhǎng)度為 4 秒，采樣頻率為 100 Hz雾家。然后我們使用一個(gè)延后量為 100（即 1 秒）的 VAR 模型與 GP CaKe 模型恢復(fù)這個(gè)脈沖響應(yīng)函數(shù)铃彰。GP CaKe 有三個(gè)重要參數(shù)，分別反映了其響應(yīng)函數(shù)的時(shí)間平滑度芯咧、時(shí)間定位和噪聲水平牙捉，我們會(huì)在下一篇文章詳細(xì)討論它們。現(xiàn)在敬飒，我直接用人工方式將這些參數(shù)設(shè)置成合理的值邪铲；在實(shí)際應(yīng)用時(shí)，我們會(huì)根據(jù)數(shù)據(jù)來估計(jì)它們的值驶拱，并通過相關(guān)應(yīng)用的背景知識(shí)來設(shè)置它們。圖 3 展示了模擬實(shí)驗(yàn)的結(jié)果晶衷±陡伲可以看到，這兩種方法都可以很好地區(qū)分當(dāng)前存在的和不存在的連接（注意圖中的縱軸是不一樣的）晌纫。對(duì)于當(dāng)前存在的連接税迷，這兩種方法都能在一定程度上恢復(fù)它的形狀；但 GP CaKe 的結(jié)果比 VAR 的結(jié)果更加平滑且更為可靠锹漱。另外箭养，這個(gè)響應(yīng)函數(shù)在 1 秒之后沒有突然的斷點(diǎn)。

圖 3：延后量為 99 的 VAR 模型與 GP CaKe 所恢復(fù)的因果脈沖響應(yīng)函數(shù)哥牍。紅線是基本真值交互毕泌，綠線是在 100 個(gè)樣本上平均后得到的恢復(fù)結(jié)果的期望喝检。灰綠色區(qū)域是 95%?置信區(qū)間撼泛。

這個(gè)模擬實(shí)驗(yàn)為 GP CaKe 在實(shí)際數(shù)據(jù)上的應(yīng)用提供了一個(gè)很好的起點(diǎn)挠说。我們看到 GP CaKe 的結(jié)果要平滑得多，也可靠得多愿题。這確實(shí)需要我們學(xué)習(xí)能確定響應(yīng)函數(shù)的平滑度损俭、定位和噪聲水平的超參數(shù)。我們的下一篇文章將回到這個(gè)主題潘酗，并還會(huì)講解該響應(yīng)函數(shù)的實(shí)際計(jì)算方法杆兵！