在MCMC之馬爾可夫鏈之中我們介紹到，給定一個概率分佈π，很難直接找到對應的馬爾可夫鏈狀態轉移矩陣P。只要解決這個問題，我們便可以找到一種通用的概率分佈取樣方法，進而用於蒙特卡羅模擬。下面我們來介紹如何找到馬爾可夫鏈所對應的狀態轉移矩陣P。

1.馬爾可夫鏈細緻平穩條件

解決平穩分佈π所對應的馬爾可夫鏈狀態轉移矩陣P之前，我們先看一下馬爾可夫鏈的細緻平穩條件。其定義為：如果非週期馬爾可夫鏈的狀態轉移矩陣P和概率分佈π(x)對於所有的i,j滿足下列方程，則概率分佈π(x)是狀態轉移矩陣P的平穩分佈。
$\pi$

( i ) P ( i , j ) = π ( j ) P ( j , i ) \pi(i)P(i,j) = \pi(j)P(j,i)
證明如下，由細緻平穩條件有
$\sum _{i=1}^{\infty}\pi(i) P(i,j) = \sum _{i=1} ^{\infty} \pi(j) P(j,i) = \pi(j) \sum _{i=1} ^{\infty}P(j,i) = \pi(j)$
將上式用矩陣表示為
$\pi P = \pi$
上式滿足馬爾可夫鏈的收斂性質，也就是說，只要我們找到可以使概率分佈π(x)滿足細緻平穩分佈的矩陣P即可。不過僅僅從細緻平穩條件還是很難找到合適的矩陣P，比如我們的目標平穩分佈使π(x)，隨機找一個馬爾可夫鏈狀態轉移矩陣Q，他是很難滿足細緻平穩條件的，即
$\pi (i) Q(i,j) \neq \pi(j) Q(j,i)$
那麼有什麼辦法可以使這個等式相等呢？

2.MCMC取樣

由於一般情況下，目標平穩分佈π(x)和某一馬爾可夫鏈狀態轉移矩陣Q不滿足細緻平穩條件，即
$\pi (i) Q(i,j) \neq \pi(j) Q(j,i)$
我們對上式進行一些變換，使細緻平穩條件成立。方法是引入一個α(i,j)，使得上式等式能夠成立，即
$\pi (i) Q(i,j) \alpha (i,j) = \pi(j) Q(j,i) \alpha(j,i)$
問題是什麼樣的α可以使上式成立？其實很簡單，只要滿足
$\alpha (i,j) = \pi (j)Q(j,i); \ \alpha(j,i)=\pi(i)Q(i,j)$
這樣，我們便找到使分佈π(x)對應的馬爾可夫鏈狀態轉移矩陣P，滿足
$P(i,j) = Q(i,j)\alpha(i,j)$
從上面可以得到，目標矩陣P可以通過任意一個馬爾可夫鏈狀態轉移矩陣Q乘以α(i,j)得到。α(i,j)我們一般稱之為接受率，取值在[0,1]之間，可以理解為一個概率值。也就是說，目標矩陣P可以通過任意一個馬爾可夫鏈狀態轉移矩陣Q以一定的接受率得到。

其實很像我們在MCMC之蒙特卡羅方法中提到的接受-拒絕取樣，那裡是以常用分佈通過一定的接受-拒絕概率得到一個非常見分佈。這裡是通過常見的馬爾可夫鏈狀態轉移矩陣Q通過一定的接受-拒絕概率得到目標轉移矩陣P，兩者解決問題的思路是相同的。下面，我們來總結下MCMC的取樣過程

• 輸入任意選定的馬爾可夫鏈狀態轉移矩陣Q，平穩分佈π(x)，設定狀態轉移次數閾值n1，需要的樣本個數n2。
• 從任意簡單概率分佈取樣得到初始值 $x_0$。
• for t=0 to n1+n2-1
  • 從條件概率分佈 $Q(x|x_t)$中取樣得到樣本 $x_*$。
  • 從均勻分佈取樣u ~ uniform[0, 1]。
  • 如果 $u < \alpha (x_t, x_*) = \pi(x_*)Q(X_*, x_t)$，則接受轉移 $x_t -> x_*$，即 $x_{t+1} = x_*$。
  • 否則不接受轉移，即 $t = max(t-1, 0)$。

樣本集 $(x_{n1},x_{n1+1},...,x_{n1+n2-1})$即為我們需要的平穩分佈樣本集。

上述過程便是MCMC取樣理論，但很難在實際應用，為什麼呢? 因為 $\alpha(x_t,x_*)$可能非常小，比如0.1，導致大部分取樣值都被拒絕轉移，取樣效率很低。可能我們取樣可上百萬次，馬爾科夫鏈還沒有收斂。實際應用中，我們可以通過M-H取樣方法進行取樣。

3.M-H取樣

M-H取樣解決了MCMC取樣接受率過低的問題，我們首先回到MCMC取樣的細緻平穩條件
$\pi (i) Q(i,j) \alpha (i,j) = \pi(j) Q(j,i) \alpha(j,i)$
取樣效率過低的原因是α(i,j)太小，比如0.1，α(j,i)為0.2，即
$\pi (i) Q(i,j) * 0.1 = \pi(j) Q(j,i) * 0.2$
如果兩邊同時擴大5倍，細緻平穩條件仍然是滿足的，即
$\pi (i) Q(i,j) * 0.5 = \pi(j) Q(j,i) * 1$

相關推薦