MCMC(三)MCMC采樣和M-H采樣
在MCMC(二)馬爾科夫鏈中我們講到給定一個概率平穩分布$\pi$, 很難直接找到對應的馬爾科夫鏈狀態轉移矩陣$P$。而只要解決這個問題，我們就可以找到一種通用的概率分布采樣方法，進而用於蒙特卡羅模擬。本篇我們就討論解決這個問題的辦法：MCMC采樣和它的易用版M-H采樣。
1. 馬爾科夫鏈的細致平穩條件
在解決從平穩分布$\pi$, 找到對應的馬爾科夫鏈狀態轉移矩陣$P$之前，我們還需要先看看馬爾科夫鏈的細致平穩條件。定義如下：
如果非周期馬爾科夫鏈的狀態轉移矩陣$P$和概率分布$\pi(x)$對於所有的$i,j$滿足：$$\pi(i)P(i,j) = \pi(j)P(j,i)$$
則稱概率分布$\pi(x)$是狀態轉移矩陣$P$的平穩分布。
證明很簡單,由細致平穩條件有：$$\sum\limits_{i=1}^{\infty}\pi(i)P(i,j) = \sum\limits_{i=1}^{\infty} \pi(j)P(j,i) = \pi(j)\sum\limits_{i=1}^{\infty} P(j,i) = \pi(j)$$
將上式用矩陣表示即為：$$\pi P = \pi$$
即滿足馬爾可夫鏈的收斂性質。也就是說，只要我們找到了可以使概率分布$\pi(x)$滿足細致平穩分布的矩陣$P$即可。這給了我們尋找從平穩分布$\pi$, 找到對應的馬爾科夫鏈狀態轉移矩陣$P$的新思路。
不過不幸的是，僅僅從細致平穩條件還是很難找到合適的矩陣$P$。比如我們的目標平穩分布是$\pi(x)$,隨機找一個馬爾科夫鏈狀態轉移矩陣$Q$,它是很難滿足細致平穩條件的，即：$$\pi(i)Q(i,j) \neq \pi(j)Q(j,i)$$
那么如何使這個等式滿足呢？下面我們來看MCMC采樣如何解決這個問題。
2. MCMC采樣
由於一般情況下，目標平穩分布$\pi(x)$和某一個馬爾科夫鏈狀態轉移矩陣$Q$不滿足細致平穩條件，即$$\pi(i)Q(i,j) \neq \pi(j)Q(j,i)$$
我們可以對上式做一個改造，使細致平穩條件成立。方法是引入一個$\alpha(i,j)$,使上式可以取等號，即：$$\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)$$
問題是什么樣的$\alpha(i,j)$可以使等式成立呢？其實很簡單，只要滿足下兩式即可：$$\alpha(i,j) = \pi(j)Q(j,i)$$$$\alpha(j,i) = \pi(i)Q(i,j)$$
這樣，我們就得到了我們的分布$\pi(x)$對應的馬爾科夫鏈狀態轉移矩陣$P$，滿足：$$P(i,j) = Q(i,j)\alpha(i,j)$$
也就是說，我們的目標矩陣$P$可以通過任意一個馬爾科夫鏈狀態轉移矩陣$Q$乘以$\alpha(i,j)$得到。$\alpha(i,j)$我們有一般稱之為接受率。取值在[0,1]之間，可以理解為一個概率值。即目標矩陣$P$可以通過任意一個馬爾科夫鏈狀態轉移矩陣$Q$以一定的接受率獲得。這個很像我們在MCMC(一)蒙特卡羅方法第4節講到的接受-拒絕采樣，那里是以一個常用分布通過一定的接受-拒絕概率得到一個非常見分布，這里是以一個常見的馬爾科夫鏈狀態轉移矩陣$Q$通過一定的接受-拒絕概率得到目標轉移矩陣$P$,兩者的解決問題思路是類似的。
好了，現在我們來總結下MCMC的采樣過程。
1）輸入我們任意選定的馬爾科夫鏈狀態轉移矩陣$Q$，平穩分布$\pi(x)$，設定狀態轉移次數閾值$n_1$，需要的樣本個數$n_2$
2）從任意簡單概率分布采樣得到初始狀態值$x_0$
3）for $t = 0$ to $n_1 +n_2-1$:
a) 從條件概率分布$Q(x|x_t)$中采樣得到樣本$x_{*}$
b) 從均勻分布采樣$u \sim uniform[0,1] $
c) 如果$u < \alpha(x_t,x_{*}) = \pi(x_{*})Q(x_{*},x_t) $, 則接受轉移$x_t \to x_{*}$，即$x_{t+1}= x_{*}$
d) 否則不接受轉移，即$x_{t+1}= x_{t}$
樣本集$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$即為我們需要的平穩分布對應的樣本集。
上面這個過程基本上就是MCMC采樣的完整采樣理論了，但是這個采樣算法還是比較難在實際中應用，為什么呢？問題在上面第三步的c步驟，接受率這兒。由於$\alpha(x_t,x_{*}) $可能非常的小，比如0.1，導致我們大部分的采樣值都被拒絕轉移，采樣效率很低。有可能我們采樣了上百萬次馬爾可夫鏈還沒有收斂，也就是上面這個$n_1$要非常非常的大，這讓人難以接受，怎么辦呢？這時就輪到我們的M-H采樣出場了。
3. M-H采樣
M-H采樣是Metropolis-Hastings采樣的簡稱，這個算法首先由Metropolis提出，被Hastings改進，因此被稱之為Metropolis-Hastings采樣或M-H采樣
M-H采樣解決了我們上一節MCMC采樣接受率過低的問題。
我們回到MCMC采樣的細致平穩條件：$$\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)$$
我們采樣效率低的原因是$\alpha(i,j)$太小了，比如為0.1，而$\alpha(j,i)$為0.2。即：$$\pi(i)Q(i,j)\times 0.1 = \pi(j)Q(j,i)\times 0.2$$
這時我們可以看到，如果兩邊同時擴大五倍，接受率提高到了0.5，但是細致平穩條件卻仍然是滿足的，即：$$\pi(i)Q(i,j)\times 0.5 = \pi(j)Q(j,i)\times 1$$
這樣我們的接受率可以做如下改進，即：$$\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$$
通過這個微小的改造，我們就得到了可以在實際應用中使用的M-H采樣算法過程如下：
1）輸入我們任意選定的馬爾科夫鏈狀態轉移矩陣$Q$，平穩分布$\pi(x)$，設定狀態轉移次數閾值$n_1$，需要的樣本個數$n_2$
2）從任意簡單概率分布采樣得到初始狀態值$x_0$
3）for $t = 0$ to $n_1 +n_2-1$:
a) 從條件概率分布$Q(x|x_t)$中采樣得到樣本$x_{*}$
b) 從均勻分布采樣$u \sim uniform[0,1] $
c) 如果$u < \alpha(x_t,x_{*}) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$, 則接受轉移$x_t \to x_{*}$，即$x_{t+1}= x_{*}$
d) 否則不接受轉移，即$x_{t+1}= x_{t}$
樣本集$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$即為我們需要的平穩分布對應的樣本集。
很多時候，我們選擇的馬爾科夫鏈狀態轉移矩陣$Q$如果是對稱的，即滿足$Q(i,j) = Q(j,i)$,這時我們的接受率可以進一步簡化為：$$\alpha(i,j) = min\{ \frac{\pi(j)}{\pi(i)},1\}$$
4. M-H采樣實例
為了更容易理解，這里給出一個M-H采樣的實例。
完整代碼參見我的github: https://github.com/ljpzzz/machinelearning/blob/master/mathematics/mcmc_3_4.ipynb
在例子里，我們的目標平穩分布是一個均值3，標准差2的正態分布，而選擇的馬爾可夫鏈狀態轉移矩陣$Q(i,j)$的條件轉移概率是以$i$為均值,方差1的正態分布在位置$j$的值。這個例子僅僅用來讓大家加深對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采樣由於接受率計算式$ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}$的存在，在高維時需要的計算時間非常的可觀，算法效率很低。同時$\alpha(i,j)$一般小於1，有時候辛苦計算出來卻被拒絕了。能不能做到不拒絕轉移呢？
2） 由於特征維度大，很多時候我們甚至很難求出目標的各特征維度聯合分布，但是可以方便求出各個特征之間的條件概率分布。這時候我們能不能只有各維度之間條件概率分布的情況下方便的采樣呢？
Gibbs采樣解決了上面兩個問題，因此在大數據時代，MCMC采樣基本是Gibbs采樣的天下，下一篇我們就來討論Gibbs采樣。
（歡迎轉載，轉載請注明出處。歡迎溝通交流： liujianping-ok@163.com）