MCMC(四)Gibbs采樣
在MCMC(三)MCMC采樣和M-H采樣中，我們講到了M-H采樣已經可以很好的解決蒙特卡羅方法需要的任意概率分布的樣本集的問題。但是M-H采樣有兩個缺點：一是需要計算接受率，在高維時計算量大。並且由於接受率的原因導致算法收斂時間變長。二是有些高維數據，特征的條件概率分布好求，但是特征的聯合分布不好求。因此需要一個好的方法來改進M-H采樣，這就是我們下面講到的Gibbs采樣。
1. 重新尋找合適的細致平穩條件
在上一篇中，我們講到了細致平穩條件：如果非周期馬爾科夫鏈的狀態轉移矩陣$P$和概率分布$\pi(x)$對於所有的$i,j$滿足：$$\pi(i)P(i,j) = \pi(j)P(j,i)$$
則稱概率分布$\pi(x)$是狀態轉移矩陣$P$的平穩分布。
在M-H采樣中我們通過引入接受率使細致平穩條件滿足。現在我們換一個思路。
從二維的數據分布開始，假設$\pi(x_1,x_2)$是一個二維聯合數據分布，觀察第一個特征維度相同的兩個點$A(x_1^{(1)},x_2^{(1)})$和$B(x_1^{(1)},x_2^{(2)})$，容易發現下面兩式成立：$$\pi(x_1^{(1)},x_2^{(1)}) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(x_1^{(1)})\pi(x_2^{(1)}|x_1^{(1)}) \pi(x_2^{(2)} | x_1^{(1)}) $$$$\pi(x_1^{(1)},x_2^{(2)}) \pi(x_2^{(1)} | x_1^{(1)}) = \pi(x_1^{(1)}) \pi(x_2^{(2)} | x_1^{(1)})\pi(x_2^{(1)}|x_1^{(1)})$$
由於兩式的右邊相等，因此我們有：$$\pi(x_1^{(1)},x_2^{(1)}) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(x_1^{(1)},x_2^{(2)}) \pi(x_2^{(1)} | x_1^{(1)}) $$
也就是：$$\pi(A) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(B) \pi(x_2^{(1)} | x_1^{(1)}) $$
觀察上式再觀察細致平穩條件的公式，我們發現在$x_1 = x_1^{(1)}$這條直線上，如果用條件概率分布$\pi(x_2| x_1^{(1)})$作為馬爾科夫鏈的狀態轉移概率，則任意兩個點之間的轉移滿足細致平穩條件！這真是一個開心的發現，同樣的道理，在在$x_2 = x_2^{(1)}$這條直線上，如果用條件概率分布$\pi(x_1| x_2^{(1)})$作為馬爾科夫鏈的狀態轉移概率，則任意兩個點之間的轉移也滿足細致平穩條件。那是因為假如有一點$C(x_1^{(2)},x_2^{(1)})$,我們可以得到：$$\pi(A) \pi(x_1^{(2)} | x_2^{(1)}) = \pi(C) \pi(x_1^{(1)} | x_2^{(1)}) $$
基於上面的發現，我們可以這樣構造分布$\pi(x_1,x_2)$的馬爾可夫鏈對應的狀態轉移矩陣$P$：$$P(A \to B) = \pi(x_2^{(B)}|x_1^{(1)})\;\; if\; x_1^{(A)} = x_1^{(B)} =x_1^{(1)}$$$$P(A \to C) = \pi(x_1^{(C)}|x_2^{(1)})\;\; if\; x_2^{(A)} = x_2^{(C)} =x_2^{(1)}$$$$P(A \to D) = 0\;\; else$$
有了上面這個狀態轉移矩陣，我們很容易驗證平面上的任意兩點$E,F$，滿足細致平穩條件：$$\pi(E)P(E \to F) = \pi(F)P(F \to E)$$
2. 二維Gibbs采樣
利用上一節找到的狀態轉移矩陣，我們就得到了二維Gibbs采樣，這個采樣需要兩個維度之間的條件概率。具體過程如下：
1）輸入平穩分布$\pi(x_1,x_2)$，設定狀態轉移次數閾值$n_1$，需要的樣本個數$n_2$
2）隨機初始化初始狀態值$x_1^{(0)}$和$x_2^{(0)}$
3）for $t = 0$ to $n_1 +n_2-1$:
a) 從條件概率分布$P(x_2|x_1^{(t)})$中采樣得到樣本$x_2^{t+1}$
b) 從條件概率分布$P(x_1|x_2^{(t+1)})$中采樣得到樣本$x_1^{t+1}$
樣本集$\{(x_1^{(n_1)}, x_2^{(n_1)}), (x_1^{(n_1+1)}, x_2^{(n_1+1)}), ..., (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)})\}$即為我們需要的平穩分布對應的樣本集。
整個采樣過程中，我們通過輪換坐標軸，采樣的過程為：$$(x_1^{(1)}, x_2^{(1)}) \to (x_1^{(1)}, x_2^{(2)}) \to (x_1^{(2)}, x_2^{(2)}) \to ... \to (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)})$$
用下圖可以很直觀的看出，采樣是在兩個坐標軸上不停的輪換的。當然，坐標軸輪換不是必須的，我們也可以每次隨機選擇一個坐標軸進行采樣。不過常用的Gibbs采樣的實現都是基於坐標軸輪換的。
3. 多維Gibbs采樣
上面的這個算法推廣到多維的時候也是成立的。比如一個n維的概率分布$\pi(x_1,x_2,...x_n)$，我們可以通過在n個坐標軸上輪換采樣，來得到新的樣本。對於輪換到的任意一個坐標軸$x_i$上的轉移，馬爾科夫鏈的狀態轉移概率為$P(x_i|x_1,x_2,...,x_{i-1},x_{i+1},...,x_n)$，即固定$n-1$個坐標軸，在某一個坐標軸上移動。
具體的算法過程如下：
1）輸入平穩分布$\pi(x_1,x_2，...,x_n)$或者對應的所有特征的條件概率分布，設定狀態轉移次數閾值$n_1$，需要的樣本個數$n_2$
2）隨機初始化初始狀態值$(x_1^{(0)},x_2^{(0)},...,x_n^{(0)}$
3）for $t = 0$ to $n_1 +n_2-1$:
a) 從條件概率分布$P(x_1|x_2^{(t)}, x_3^{(t)},...,x_n^{(t)})$中采樣得到樣本$x_1^{t+1}$
b) 從條件概率分布$P(x_2|x_1^{(t+1)}, x_3^{(t)}, x_4^{(t)},...,x_n^{(t)})$中采樣得到樣本$x_2^{t+1}$
c)...
d) 從條件概率分布$P(x_j|x_1^{(t+1)}, x_2^{(t+1)},..., x_{j-1}^{(t+1)},x_{j+1}^{(t)}...,x_n^{(t)})$中采樣得到樣本$x_j^{t+1}$
e)...
f) 從條件概率分布$P(x_n|x_1^{(t+1)}, x_2^{(t+1)},...,x_{n-1}^{(t+1)})$中采樣得到樣本$x_n^{t+1}$
樣本集$\{(x_1^{(n_1)}, x_2^{(n_1)},..., x_n^{(n_1)}), ..., (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)},...,x_n^{(n_1+n_2-1)})\}$即為我們需要的平穩分布對應的樣本集。
整個采樣過程和Lasso回歸的坐標軸下降法算法非常類似，只不過Lasso回歸是固定$n-1$個特征，對某一個特征求極值。而Gibbs采樣是固定$n-1$個特征在某一個特征采樣。
同樣的，輪換坐標軸不是必須的，我們可以隨機選擇某一個坐標軸進行狀態轉移，只不過常用的Gibbs采樣的實現都是基於坐標軸輪換的。
4. 二維Gibbs采樣實例
這里給出一個Gibbs采樣的例子。完整代碼參見我的github: https://github.com/ljpzzz/machinelearning/blob/master/mathematics/mcmc_3_4.ipynb
假設我們要采樣的是一個二維正態分布$Norm(\mu,\Sigma)$,其中：$$\mu = (\mu_1,\mu_2) = (5,-1)$$$$\Sigma = \left( \begin{array}{ccc}
\sigma_1^2&\rho\sigma_1\sigma_2 \\
\rho\sigma_1\sigma_2 &\sigma_2^2 \end{array} \right) = \left( \begin{array}{ccc}
1&1 \\
1&4 \end{array} \right)$$
而采樣過程中的需要的狀態轉移條件分布為：$$P(x_1|x_2) = Norm\left ( \mu _1+\rho \sigma_1/\sigma_2 \left ( x _2-\mu _2 \right ), (1-\rho ^2)\sigma_1^2 \right )$$$$P(x_2|x_1) = Norm\left ( \mu _2+\rho \sigma_2/\sigma_1 \left ( x _1-\mu _1 \right ), (1-\rho ^2)\sigma_2^2 \right )$$
具體的代碼如下：
from mpl_toolkits.mplot3d import Axes3D from scipy.stats import multivariate_normal samplesource = multivariate_normal(mean=[5,-1], cov=[[1,1],[1,4]]) def p_ygivenx(x, m1, m2, s1, s2): return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt((1 - rho ** 2) * (s2**2)))) def p_xgiveny(y, m1, m2, s1, s2): return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt((1 - rho ** 2) * (s1**2)))) N = 5000 K = 20 x_res = [] y_res = [] z_res = [] m1 = 5 m2 = -1 s1 = 1 s2 = 2 rho = 0.5 y = m2 for i in xrange(N): for j in xrange(K): x = p_xgiveny(y, m1, m2, s1, s2) y = p_ygivenx(x, m1, m2, s1, s2) z = samplesource.pdf([x,y]) x_res.append(x) y_res.append(y) z_res.append(z) num_bins = 50 plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5) plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5) plt.title('Histogram') plt.show()
輸出的兩個特征各自的分布如下：
然后我們看看樣本集生成的二維正態分布，代碼如下：
fig = plt.figure() ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20) ax.scatter(x_res, y_res, z_res,marker='o') plt.show()
輸出的正態分布圖如下：
5. Gibbs采樣小結
由於Gibbs采樣在高維特征時的優勢，目前我們通常意義上的MCMC采樣都是用的Gibbs采樣。當然Gibbs采樣是從M-H采樣的基礎上的進化而來的，同時Gibbs采樣要求數據至少有兩個維度，一維概率分布的采樣是沒法用Gibbs采樣的,這時M-H采樣仍然成立。
有了Gibbs采樣來獲取概率分布的樣本集，有了蒙特卡羅方法來用樣本集模擬求和，他們一起就奠定了MCMC算法在大數據時代高維數據模擬求和時的作用。MCMC系列就在這里結束吧。
（歡迎轉載，轉載請注明出處。歡迎溝通交流： liujianping-ok@163.com）