MCMC(三)MCMC採樣和M-H採樣
1. 馬爾科夫鏈的細緻平穩條件
在解決從平穩分佈ππ, 找到對應的馬爾科夫鏈狀態轉移矩陣PP之前，我們還需要先看看馬爾科夫鏈的細緻平穩條件。定義如下：
如果非週期馬爾科夫鏈的狀態轉移矩陣PP和概率分佈π(x)π(x)對於所有的i,ji,j滿足：
則稱概率分佈π(x)π(x)是狀態轉移矩陣PP的平穩分佈。
證明很簡單,由細緻平穩條件有：
將上式用矩陣表示即爲：
即滿足馬爾可夫鏈的收斂性質。也就是說，只要我們找到了可以使概率分佈π(x)π(x)滿足細緻平穩分佈的矩陣PP即可。這給了我們尋找從平穩分佈ππ, 找到對應的馬爾科夫鏈狀態轉移矩陣PP的新思路。
不過不幸的是，僅僅從細緻平穩條件還是很難找到合適的矩陣PP。比如我們的目標平穩分佈是π(x)π(x),隨機找一個馬爾科夫鏈狀態轉移矩陣QQ,它是很難滿足細緻平穩條件的，即：
那麼如何使這個等式滿足呢？下面我們來看MCMC採樣如何解決這個問題。
2. MCMC採樣
由於一般情況下，目標平穩分佈π(x)π(x)和某一個馬爾科夫鏈狀態轉移矩陣QQ不滿足細緻平穩條件，即
我們可以對上式做一個改造，使細緻平穩條件成立。方法是引入一個α(i,j)α(i,j),使上式可以取等號，即：
問題是什麼樣的α(i,j)α(i,j)可以使等式成立呢？其實很簡單，只要滿足下兩式即可：
這樣，我們就得到了我們的分佈π(x)π(x)對應的馬爾科夫鏈狀態轉移矩陣PP，滿足：
也就是說，我們的目標矩陣PP可以通過任意一個馬爾科夫鏈狀態轉移矩陣QQ乘以α(i,j)α(i,j)得到。α(i,j)α(i,j)我們有一般稱之爲接受率。取值在[0,1]之間，可以理解爲一個概率值。即目標矩陣PP可以通過任意一個馬爾科夫鏈狀態轉移矩陣QQ以一定的接受率獲得。這個很像我們在MCMC(一)蒙特卡羅方法第4節講到的接受-拒絕採樣，那裏是以一個常用分佈通過一定的接受-拒絕概率得到一個非常見分佈，這裏是以一個常見的馬爾科夫鏈狀態轉移矩陣QQ通過一定的接受-拒絕概率得到目標轉移矩陣PP,兩者的解決問題思路是類似的。
好了，現在我們來總結下MCMC的採樣過程。
1）輸入我們任意選定的馬爾科夫鏈狀態轉移矩陣QQ，平穩分佈π(x)π(x)，設定狀態轉移次數閾值n1n1，需要的樣本個數n2n2
2）從任意簡單概率分佈採樣得到初始狀態值x0x0
3）for t=0t=0 to n1+n2−1n1+n2−1:
a) 從條件概率分佈Q(x|xt)Q(x|xt)中採樣得到樣本x∗x∗
b) 從均勻分佈採樣u∼uniform[0,1]u∼uniform[0,1]
c) 如果u<α(xt,x∗)=π(x∗)Q(x∗,xt)u<α(xt,x∗)=π(x∗)Q(x∗,xt), 則接受轉移xt→x∗xt→x∗，即xt+1=x∗xt+1=x∗
d) 否則不接受轉移，t=max(t−1,0)t=max(t−1,0)
樣本集(xn1,xn1+1,...,xn1+n2−1)(xn1,xn1+1,...,xn1+n2−1)即爲我們需要的平穩分佈對應的樣本集。
上面這個過程基本上就是MCMC採樣的完整採樣理論了，但是這個採樣算法還是比較難在實際中應用，爲什麼呢？問題在上面第三步的c步驟，接受率這兒。由於α(xt,x∗)α(xt,x∗)可能非常的小，比如0.1，導致我們大部分的採樣值都被拒絕轉移，採樣效率很低。有可能我們採樣了上百萬次馬爾可夫鏈還沒有收斂，也就是上面這個n1n1要非常非常的大，這讓人難以接受，怎麼辦呢？這時就輪到我們的M-H採樣出場了。
3. M-H採樣
M-H採樣是Metropolis-Hastings採樣的簡稱，這個算法首先由Metropolis提出，被Hastings改進，因此被稱之爲Metropolis-Hastings採樣或M-H採樣
M-H採樣解決了我們上一節MCMC採樣接受率過低的問題。
我們回到MCMC採樣的細緻平穩條件：
我們採樣效率低的原因是α(i,j)α(i,j)太小了，比如爲0.1，而α(j,i)α(j,i)爲0.2。即：
這時我們可以看到，如果兩邊同時擴大五倍，接受率提高到了0.5，但是細緻平穩條件卻仍然是滿足的，即：
這樣我們的接受率可以做如下改進，即：
通過這個微小的改造，我們就得到了可以在實際應用中使用的M-H採樣算法過程如下：
1）輸入我們任意選定的馬爾科夫鏈狀態轉移矩陣QQ，平穩分佈π(x)π(x)，設定狀態轉移次數閾值n1n1，需要的樣本個數n2n2
2）從任意簡單概率分佈採樣得到初始狀態值x0x0
3）for t=0t=0 to n1+n2−1n1+n2−1:
a) 從條件概率分佈Q(x|xt)Q(x|xt)中採樣得到樣本x∗x∗
b) 從均勻分佈採樣u∼uniform[0,1]u∼uniform[0,1]
c) 如果u<α(xt,x∗)=min{π(j)Q(j,i)π(i)Q(i,j),1}u<α(xt,x∗)=min{π(j)Q(j,i)π(i)Q(i,j),1}, 則接受轉移xt→x∗xt→x∗，即xt+1=x∗xt+1=x∗
d) 否則不接受轉移，t=max(t−1,0)t=max(t−1,0)
樣本集(xn1,xn1+1,...,xn1+n2−1)(xn1,xn1+1,...,xn1+n2−1)即爲我們需要的平穩分佈對應的樣本集。
很多時候，我們選擇的馬爾科夫鏈狀態轉移矩陣QQ如果是對稱的，即滿足Q(i,j)=Q(j,i)Q(i,j)=Q(j,i),這時我們的接受率可以進一步簡化爲：
4. M-H採樣實例
爲了更容易理解，這裏給出一個M-H採樣的實例。
在例子裏，我們的目標平穩分佈是一個均值3，標準差2的正態分佈，而選擇的馬爾可夫鏈狀態轉移矩陣Q(i,j)Q(i,j)的條件轉移概率是以ii爲均值,方差1的正態分佈在位置jj的值。這個例子僅僅用來讓大家加深對M-H採樣過程的理解。畢竟一個普通的一維正態分佈用不着去用M-H採樣來獲得樣本。
代碼如下：
import random import math from scipy.stats import norm import matplotlib.pyplot as plt %matplotlib inline def norm_dist_prob(theta): y = norm.pdf(theta, loc=3, scale=2) return y T = 5000 pi = [0 for i in range(T)] sigma = 1 t = 0 while t < T-1: t = t + 1 pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None) alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1]))) u = random.uniform(0, 1) if u < alpha: pi[t] = pi_star[0] else: pi[t] = pi[t - 1] plt.scatter(pi, norm.pdf(pi, loc=3, scale=2)) num_bins = 50 plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7) plt.show()
輸出的圖中可以看到採樣值的分佈與真實的分佈之間的關係如下，採樣集還是比較擬合對應分佈的。
5. M-H採樣總結
M-H採樣完整解決了使用蒙特卡羅方法需要的任意概率分佈樣本集的問題，因此在實際生產環境得到了廣泛的應用。
但是在大數據時代，M-H採樣面臨着兩大難題：
1） 我們的數據特徵非常的多，M-H採樣由於接受率計算式π(j)Q(j,i)π(i)Q(i,j)π(j)Q(j,i)π(i)Q(i,j)的存在，在高維時需要的計算時間非常的可觀，算法效率很低。同時α(i,j)α(i,j)一般小於1，有時候辛苦計算出來卻被拒絕了。能不能做到不拒絕轉移呢？
2） 由於特徵維度大，很多時候我們甚至很難求出目標的各特徵維度聯合分佈，但是可以方便求出各個特徵之間的條件概率分佈。這時候我們能不能只有各維度之間條件概率分佈的情況下方便的採樣呢？
Gibbs採樣解決了上面兩個問題，因此在大數據時代，MCMC採樣基本是Gibbs採樣的天下，下一篇我們就來討論Gibbs採樣。
本文轉自劉建平Pinard博客園博客，原文鏈接：http://www.cnblogs.com/pinard/p/6638955.html，如需轉載請自行聯繫原作者