MCMC採樣和M-H採樣

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

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

解決平穩分佈π所對應的馬爾可夫鏈狀態轉移矩陣P之前，我們先看一下馬爾可夫鏈的細緻平穩條件。其定義爲：如果非週期馬爾可夫鏈的狀態轉移矩陣P和概率分佈π(x)對於所有的i,j滿足下列方程，則概率分佈π(x)是狀態轉移矩陣P的平穩分佈。
$\pi(i)P(i,j) = \pi(j)P(j,i)$π(i)P(i,j)=π(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)$i=1∑∞​π(i)P(i,j)=i=1∑∞​π(j)P(j,i)=π(j)i=1∑∞​P(j,i)=π(j)
將上式用矩陣表示爲
$\pi P = \pi$πP=π
上式滿足馬爾可夫鏈的收斂性質，也就是說，只要我們找到可以使概率分佈π(x)滿足細緻平穩分佈的矩陣P即可。不過僅僅從細緻平穩條件還是很難找到合適的矩陣P，比如我們的目標平穩分佈使π(x)，隨機找一個馬爾可夫鏈狀態轉移矩陣Q，他是很難滿足細緻平穩條件的，即
$\pi (i) Q(i,j) \neq \pi(j) Q(j,i)$π(i)Q(i,j)̸​=π(j)Q(j,i)
那麼有什麼辦法可以使這個等式相等呢？

2.MCMC採樣

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

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

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

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

上述過程便是MCMC採樣理論，但很難在實際應用，爲什麼呢? 因爲 $\alpha(x_t,x_*)$α(xt​,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)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)
採樣效率過低的原因是α(i,j)太小，比如0.1，α(j,i)爲0.2，即
$\pi (i) Q(i,j) * 0.1 = \pi(j) Q(j,i) * 0.2$π(i)Q(i,j)∗0.1=π(j)Q(j,i)∗0.2
如果兩邊同時擴大5倍，細緻平穩條件仍然是滿足的，即
$\pi (i) Q(i,j) * 0.5 = \pi(j) Q(j,i) * 1$π(i)Q(i,j)∗0.5=π(j)Q(j,i)∗1
這樣我們可以對接受率做如下改進，即
$\alpha (i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1 \}$α(i,j)=min{π(i)Q(i,j)π(j)Q(j,i)​,1}
通過上述的轉換，我們便可在實際應用中使用M-H算法進行採樣，M-H採樣算法過程如下所示

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

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

很多時候，我們選擇的馬爾科夫鏈狀態轉移矩陣Q如果是對稱的，即滿足Q(i,j)=Q(j,i)，這時我們的接受率可以進一步簡化爲
$\alpha (i,j) = min(\frac{\pi(j)}{\pi(i)},1)$α(i,j)=min(π(i)π(j)​,1)

4.M-H採樣總結

M-H採樣解決了使用蒙特卡羅方法需要的任意概率分佈樣本集的問題，因此在實際生產環境中得到廣泛應用。但在大數據情況下，M-H面臨如下問題

• **數據特徵非常多：**因爲M-H採樣由於接受率 $\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}$π(i)Q(i,j)π(j)Q(j,i)​的存在，在高維計算時需要很長的計算時間，算法效率很低。同時α(i,j)一般小於1，有時候辛苦計算出來的結果卻被拒絕，能不能做到不拒絕轉移呢？
• **特徵維度比較大：**很多時候我們很難求出目標的各特徵維度聯合分佈，但是可以方便求出各個特徵之間的條件概率分佈。這時能不能只用各維度之間的條件概率分佈去方便的採樣呢?

5.推廣

更多內容請關注公衆號謂之小一，若有疑問可在公衆號後臺提問，隨時回答，歡迎關注，內容轉載請註明出處。