Dynamic Causal Modeling:確定性因果模型(下)

分類:實用技巧 時間:2017-09-24

對於(上)的小結:

  • 從心理學刺激輸入到bold信號建立起了一套非線性的數學模型,這個數學模型的參數可以分為兩類: A,Bj和C是描述神經元活動屬性的參數,統稱為θc, τi、κi、Υi、α的ρ是描述血管動力學屬性的參數,統稱為θh。 θc和θh統稱為θ。設y為fMRI儀器讀取的bold信號,u為刺激輸入,那麽得到描述刺激輸入和bold信號的非線性方程:
  • 上式中,我們已經知道了u和y,需要估計參數θ,即一系列包含大腦連接變化的參數。DCM主要使用貝葉斯參數估計的方法來估計θ。

第一部分 貝葉斯參數估計:

頻率理論:樣本x=(x1,…,xn)的產生只需要一個步驟,即人們從總體中抽取了一系列樣本。我們可以直接從這些樣本值估計出總體的參數。

貝葉斯理論:樣本x=(x1,…,xn)的產生需要兩個步驟:

  • 第一步,從總體的參數分布(先驗分布)中,由“上天”抽取了一個參數θ’;
  • 第二步,這個參數確定後的總體的分布中,抽取了一系列樣本。

因此,我們不但需要從樣本值做參數的“無偏估計”,還要根據參數本身的分布規律(先驗分布)讓無偏估計值“偏”一下(在此不評論,也不做各種解釋,形象直觀地描述一下)。

普通概率論課本中貝葉斯公式(這個學過概率論的同學都已經熟悉了):

對等式右側分母(即m(A):A的邊緣密度函數)僅僅起到正則化的作用,因此可以得到下面的正比關系,註意到這裏面已經將A、B換成了y和θ,即當樣本值為y時,參數θ的分布規律可由下式來確定。

在貝葉斯統計理論中,p(θ|y)是已知的樣本量(比如說bold信號)決定的參數(比如說大腦連接強度)的分布規律。被稱為後驗概率。

p(θ)是人們對於參數分布的先驗假設。即先驗分布。

p(y|θ)是人們對於樣本分布的似然函數。概率論通常用這個分布估計參數。而貝葉斯理論認為還要添加先驗概率的影響。

上面這個式子可以簡單概括為:後驗概率就是極大似然(無偏估計)和先驗概率的折中方案。

第二部分 參數先驗分布的確定

DCM模型中,待估參數有8個,其中,血管動力學參數的先驗分布已經由Friston等人在2002年的研究中確定,具體值見附錄。

而神經活動的參數A(靜息態的腦區間連接)、B(刺激引起的腦區間連接變化)、C(輸入引起的目標腦區神經活動變化)的先驗分布,則要通過構建一個穩定的網絡結構來確定。

(1)穩定系統和非穩定系統

  • 穩定性是系統的一種特征。涉及到在時間演化過程中,兩個在相平面上相鄰無限近的點會不會必然分離。如下圖就是一個極端不穩定的系統(混沌系統),系統初始點非常鄰近,但隨著時間進展出現了分叉和混沌現象。DCM方法通過對系統的穩定性進行先驗性假設,來形成對A、B、C三個參數的先驗分布。

很明顯,神經活動的狀態z不會隨著時間而指數發散,而是會回到一個穩定的狀態,這就意味著固有連接矩陣A的最大實特征根為負值,即矩陣A的主李雅普諾夫指數應為負值(這是非線性動力學判定系統穩定性的方法之一,如有疑問請閱讀附錄)。

(2)通過相平面認識系統穩定性:

先將下圖中x和x’想象成鐘擺的重力勢能和動能。實際上相平面就是勢能和動能的坐標圖。隨著時間演進,系統中的勢能和動能會互相轉化。造成了點在相平面上運動。

當x’=0時,系統陷入不動點(奇點,定點),即陷入停滯狀態。一般通過李雅普諾夫指數(方法)判定系統解的穩定性 ,在此定義下,李雅普諾夫穩定性指的是系統在想平面上的運動軌跡受到微擾後,是否還會回到原來的軌道上。人們通過李雅普諾夫指數來檢驗和描述系統的穩定性。(確切定義和擴展閱讀看附錄) 。

(3)通過系統穩定性假設得到參數先驗分布

回到DCM模型中來,我們顯然可以假設:人的靜息態大腦連接(A)這一個系統不會陷入到劇烈變動狀態,刺激輸入後,大腦連接的變動(B)也不會無限擴散,外界刺激對於大腦激活(C)也會隨著時間進展慢慢恢復到靜息態。

如果我們先對A、B、C三個參數進行標準化,那麽它們標準化之後的絕對零點,同時也是高維相平面上的絕對零點。因此,只要確保A、B、C三個矩陣的主李雅普諾夫指數最小,就可以保證系統是穩定的。這是我們對於這三個參數的先驗假設。

對於連接矩陣A,設一共有l個腦區,從i腦區到j腦區的連接度為aij, ζ為它們的平方和,當它們都相等時,得到主李雅普諾夫λa指數最大,這兩個值為:

根據正態分布的特點,可以認為這兩個式子(在絕對零點附近)保持局部線性,因為我們要保證λa<0,這樣我們就建立起高斯分布重要的參數平方和和系統穩定性之間的關系: ζ<l/l-1。

設aij符合一個均值為0,方差為va的正態分布,則aij的平方和符合卡方分布,在這個分布假設下,對於參數進行估計後,就可以得到A矩陣元素分布的超參數——方差:

這樣,我們就對於A的分布有了一些先驗性的看法,對A標準化後,它是一個均值為0,方差為v的(取0- va的無假設分布)的正態分布。

B矩陣的元素也可以用同樣的方法得到先驗分布。對於C矩陣,Friston等人假設的是它的元素符合均值為0,標準差為1的正態分布。

第三部分 數值求解

現在已經得到了DCM的基本模型,從外界刺激輸入u到bold信號,建立起了一套非線性關系。還有了關於所有待估計參數的先驗分布。下面將根據它們的先驗分布,采用貝葉斯參數估計的方法將這些參數的後驗分布求出來。

現在有兩個困難:

  • 1、DCM模型中,輸入u和bold信號y之間的關系是非線性的,那麽比如基於最小二乘回歸模型等要求基本線性假設的參數估計方法就不能使用了。
  • 2、如果用直接解極大似然法構造的積分方程來估計這些參數,那麽又會遇到最大似然估計容易遇到的問題:很難找到通用法則來計算結果。

Friston等人首先利局部線性近似的方法進行構造參數估計的通用法則。

局部線性近似:指的是穩定系統在某點附近處運動時,運動的速率近似於平衡點的運動速率。也因此,我們可以在參數θ分布的中心(眾數、均值)附近以直代曲,認為在這一點附近是線性變化。

根據局部線性近似的思想,Friston認為各個參數都是圍繞著均值上下微動,斜率和均值點附近的一致。設在i腦區參數均值為

根據貝葉斯參數估計方法,假設後驗分布、似然函數都是正態分布:

Firston在2002年指出,在層級觀察模型中(Hierarchical Linear Observational

Model)兩個層級(region to region)的壓縮到一個非層級方程中時(沒有找到專門的漢語翻譯,可能我說的不太專業,歡迎有統計背景的同學指出錯誤),可以得到以下關系:

?將上式代換到前面的先驗和最大似然等式中去,又根據貝葉斯參數估計方法,得到:

上述推導中,假設測量誤差ε的協方差Cε已知,但其實是未知的,根據Neal和Heaton(1998)的研究,對於協方差C ε的估計,可以找到一個適用於所有樣本類型的超參數估計方法,運用樣本值估計測量誤差的分布:

到目前為止,已經找到了根據刺激輸入u、bold信號、第一個腦區的參數θ(1)去估計其它所有腦區的參數θ(i)的方法。

DCM模型假定,第一個腦區的參數就是前面提到的先驗分布中確定的參數均值。

第三部分 算法和步驟

  • 1、選定ROI腦區,確定連接和連接方向,以及各腦區是否受到外界刺激的直接影響。
  • 2、將模型中第一個腦區的參數值θ和λ設置為先驗均值。
  • 3、采用EM算法,推算下一個腦區的參數θ和λ的後驗分布的方差和均值,直到將所有腦區的參數都計算完畢。
  • 4、顯著性檢驗:估計得到了如腦區連接變化矩陣Bj的方差和均值(隱式),進行拉普拉斯近似,得到後驗分布,給定某個閾值如α,直接看α是否落在後驗分布的拒絕區域。
  • 5、驗證不同模型,通過計算模型的貝葉斯因子來比較哪個模型比較好。

附錄:

1、Friston在spm中規定的血管動力學參數的先驗分布:

κ:mean=0.65/s sd=0.015

γ: mean=0.41/s sd=0.002

τ:mean=0.98/s sd=0.0568

α: mean=0.32 sd=0.0015

ρ:mean=0.34 sd=0.0024

2、李雅普諾夫穩定性和判定(因為是在PPT整理的,直接傳圖)

正如以前說的,DCM的原理比較復雜,有計算神經科學基礎的同學如想在本文之外做進一步閱讀,可以看胡德文老師等寫的《腦磁共振影像數據的時空分析》中的專門介紹章節(其中沒有介紹層級模型,而且有些推導過程太跳躍,不適合沒有相關背景的同學閱讀,可結合本文看)。其中,拉普拉斯近似、EM算法等算法內容大家都可以百度到,就不作贅述。

公式是我從論文摳圖下來或者自己在Office中輸入摳圖的。大小不一。實在沒有時間再整理一遍,也是無可奈何之舉。

如有錯誤之處希望大家指教!


Tags: 參數 貝葉 樣本 分布 估計 統稱

文章來源:


ads
ads

相關文章
ads

相關文章

ad