使用高斯過程的因果推理：GP CaKe 的基本思路

選自mindcodec

機器之心編譯

參與：Panda

荷蘭拉德堡德大學的研究者的一篇 NIPS 2017 論文基於對大腦中因果交互的研究，提出了一種用於因果推理的通用方法 GP CaKe。日前，相關研究者開始通過一系列文章講解這種融合了向量自回歸模型與動態因果建模兩者之長的新方法的思路及其應用。本文為該系列的第一篇，闡釋了 GP CaKe 的設計思路。

項目地址：https://github.com/LucaAmbrogioni/GP-CaKe-project

原論文：http://papers.nips.cc/paper/6696-gp-cake-effective-brain-connectivity-with-causal-kernels.pdf

我們最近開發出了一種用於時間序列數據中因果推理的全新方法 [Ambrogioni et al., 2017]。我們稱之為「GP CaKe」，即具有因果核的高斯過程（Gaussian Processes with Causal Kernels）。這種方法不僅縮寫詞（有「蛋糕」的意思）很美味，而且將向量自回歸模型（VAR）那引人注目的特性與動態因果建模（DCM）優雅地結合到了一起。是的，確實是一舉兩得！

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

背景：對多變數時間序列的分析

這項研究的研究背景是帶有一個時間維度的複雜系統。在我們的案例中可能是對多個大腦區域的激活狀態的連續測量，比如通過 EEG、MEG 或 fMRI [2] 等方式，但也可以是連續的股票交易列表、天氣現象以及蛋白質濃度變化等等。

在統計學和機器學習領域，建模這種複雜系統的時間序列的方法主要有兩種：向量自回歸（VAR）[Lütkepohl, 2005] 和動態系統理論（DST），後者通常是通過（隨機）微分方程（SDE）或差分方程 [Abraham & Shaw, 1983] 實現的。我們將簡要介紹這兩種方法，以激發你的興趣以及說明 GP CaKe 與它們的關聯。如果你已經很了解這兩種方法了，可以直接跳過這部分，查看後文有關 GP CaKe 的內容。

向量自回歸

向量自回歸（VAR）最基本的形式是預測一個特定變數 x_j(t) 在時間 t 的值，該變數是另一些變數 x_i 的（隨機）函數，它們的關係如下 [3]：

這個等式的含義如下：

變數 x_j(t) 的信號取決於該變數從所有其它變數獲得的輸入。這種依賴關係的強度由自回歸係數 a_ij(τ) 確定。參數 τ 是信號 x_j(t) 和 x_i(t) 之間的延後量（lag）。綜合起來，這表示一個變數對另一個變數的影響可以在 τ=0 時為零（舉個例子），然後這種因果影響會緩慢增大（即 a_ij(τ) 更大），只有當延後量變得很大時才會再次衰減——這意味著在遙遠過去發生的某些事情現在已經不再重要了。如果我們繪製出這些係數與該延後量的函數關係，我們就會得到所謂的脈衝響應函數（impulse response function/IRF）。最後，w(t) 描述的是驅動系統的隨機「創新量（innovations）」或「衝擊量（shocks）」。它們能夠反映 x_j(t) 的內部動態。比如說，我們國家的天氣會受周邊鄰國天氣狀況的影響（即它們對我們的氣候有因果影響），但也會受到我國內部情況的影響。如果 a_ij(τ)>0，我們就可以說 x_i 對 x_j 有因果影響（這一思想的實際實現將需要某些顯著性測試）。這意味著因果關係的時間概念：一個變數的過去能為我們提供有關另一個變數的未來的信息。這種看待因果關係的角度也被稱為 Wiener-Granger 因果關係，有時也被簡稱為 Granger 因果關係 [Bressler & Seth, 2011]。通過觀察 IRF，我們可以明確地看到 Granger 因果交互的時間形狀：

圖 1：在三大金融指數上的 VAR 分析示例。左圖是每個變數的時間序列，右圖是最大為 10 個月的延後量的脈衝響應函數（IRF）。注意其中也包含了自響應。

動態系統理論

顧名思義，動態系統理論（DST）也就是建模了系統的動態。以經典的 Ornstein-Uhlenbeck 過程為例，如下所示：

這描述了一個隨機遊走過程——隨著時間推移，會逐漸趨近其平均值 μ。

圖 2：五個全都（漸進地）回到同一平均值 μ=0.8 的 Ornstein-Uhlenbeck 過程，但它們各自的雜訊水平和回到均值的速度都各不相同。

DST 在動態因果建模（DCM）方面也有應用 [Friston, 2009]。儘管 DCM 的大多數實現都包含一個專用於 fMRI 的前向模型，將 DCM 的應用限制在了神經成像研究領域，但其核心有一個通用型的微分方程系統：

注意，其中 X=(x_1, ..., x_n)。此外，A 是一個包含 x 中變數之間固定交互的矩陣。其作用類似於 VAR 模型的自回歸係數，但 DCM 中沒有建模延後量。實際上，其動態所受的影響是即時的。另外的 B 和 C 項表示（節點特定的）外源輸入 μ，我們這裡不會對此詳做討論，因為 GP CaKe [4] 中（尚）沒有與這些項類同的地方。

連續和動態的向量自回歸：GP CaKe

前提內容已經足夠了，現在該進入正題了 [5]。VAR 模型的難點是在實際操作中，我們沒有足夠多的觀察來可靠地估計自回歸係數。由此造成我們的脈衝響應函數充滿雜訊且難以解讀。此外，VAR 模型只能粗略地描述系統的動態。高階交互會被完全忽視。DCM 確實能更廣泛地考慮動態，但是它卻不能建模一個變數的變化和另一個變數的動態的變化之間的延遲情況，這是有問題的。DCM 的某些變體確實包含一個延後項，但會將其設置為一個常量項，而不是我們估計交互係數的值的一個區間。你可能也已經猜到了：GP CaKe 實際上將延後的交互與動態系統結合了起來。下面來一窺究竟。

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

其中，D_j 是微分運算元（即其描述了直到第 p 個導數的動態），w_j(t) 仍然是指「創新量（innovation）」或「衝擊量（shock）」項，關鍵的 C_j(t) 是來自其它變數 i≠j 的因果效應的總和。

其中，C_i(t) 到 C_j(t) 是因果脈衝響應函數（CIRF），描述了從 i 到 j 的因果交互。從中可以看到，C_j(t) 是時間序列的求和（在所有輸入變數上），並且這些時間序列與它們的脈衝響應函數進行了卷積。這個定義完全類似於 VAR 模型等式右側的第一項，但卻是連續的，而非離散的。然而，GP CaKe 並不簡單地是 VAR 的連續式變體。微分運算元 D_j 看似沒啥作用，但實際非常關鍵。它描述的是一個變數的內部動態，且無論它從其它變數那裡得到了什麼輸入——而且我們尚未描述這些動態是什麼！存在一些（實際上數量很有限）可能的選項，比如，這些動態可以是我們之前見過的簡單的 Ornstein-Uhlenbeck 隨機遊走，即一個振蕩過程。不管是在什麼案例中，都要記住這樣一個重點：GP CaKe 假設來自其它變數的輸入會通過因果脈衝響應函數影響動態 D_jx_j(t)，而不只會影響 x_j(t) 本身！

讓我們實現它

在下一篇文章中，我將解釋如何計算因果脈衝響應函數，這與我們之前的一篇有關高斯過程回歸中傅立葉變換的文章有關：https://www.mindcodec.com/the-fourier-transform-through-the-lens-of-gaussian-process-regression/。目前，我們就假設我們已經有能幫助我們完成這項任務的工具了（實際上我們確實有，參見對應的 GitHub 項目）。我們這裡只是為了進行一點演示，不會涉及太多細節。我們將使用一個已知的脈衝響應函數來生成某些數據，然後會嘗試使用 VAR 和 GP CaKe 來恢復它。注意，目前已經存在一些比標準的非正則化的 VAR 更好的實現，但這裡只是用於解釋說明，普通的 VAR 就足夠了。

我們從兩個變數 x_A 和 x_B 開始，它們具有如下的因果關係：

其中 τ 仍然是兩個變數之間的時間延後量，而 s 則是脈衝響應的長度範圍（這個函數的形狀如下圖中的紅線所示）。為了我們當前的目的，這是一個任意的變數，我們只是隨便選取了某個值。此外，我們假設這兩個變數的內部動態都是 Ornstein-Uhlenbeck 過程，因此

其中，α 是該過程的弛豫係數，指示了該時間序列回到其均值（零）的速度。我們為該動態系統生成了 100 個樣本，總長度為 4 秒，採樣頻率為 100 Hz。然後我們使用一個延後量為 100（即 1 秒）的 VAR 模型與 GP CaKe 模型恢復這個脈衝響應函數。GP CaKe 有三個重要參數，分別反映了其響應函數的時間平滑度、時間定位和雜訊水平，我們會在下一篇文章詳細討論它們。現在，我直接用人工方式將這些參數設置成合理的值；在實際應用時，我們會根據數據來估計它們的值，並通過相關應用的背景知識來設置它們。圖 3 展示了模擬實驗的結果。可以看到，這兩種方法都可以很好地區分當前存在的和不存在的連接（注意圖中的縱軸是不一樣的）。對於當前存在的連接，這兩種方法都能在一定程度上恢復它的形狀；但 GP CaKe 的結果比 VAR 的結果更加平滑且更為可靠。另外，這個響應函數在 1 秒之後沒有突然的斷點。

圖 3：延後量為 99 的 VAR 模型與 GP CaKe 所恢復的因果脈衝響應函數。紅線是基本真值交互，綠線是在 100 個樣本上平均後得到的恢復結果的期望。灰綠色區域是 95% 置信區間。

這個模擬實驗為 GP CaKe 在實際數據上的應用提供了一個很好的起點。我們看到 GP CaKe 的結果要平滑得多，也可靠得多。這確實需要我們學習能確定響應函數的平滑度、定位和雜訊水平的超參數。我們的下一篇文章將回到這個主題，並還會講解該響應函數的實際計算方法！

本文為機器之心編譯，轉載請聯繫本公眾號獲得授權。

------------------------------------------------

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！