在MCMC(三)MCMC取樣和M-H取樣中，我們講到了M-H取樣已經可以很好的解決蒙特卡羅方法需要的任意概率分佈的樣本集的問題。但是M-H取樣有兩個缺點：一是需要計算接受率，在高維時計算量大。並且由於接受率的原因導致演算法收斂時間變長。二是有些高維資料，特徵的條件概率分佈好求，但是特徵的聯合分佈不好求。因此需要一個好的方法來改進M-H取樣，這就是我們下面講到的Gibbs取樣。

1. 重新尋找合適的細緻平穩條件

　　　　在上一篇中，我們講到了細緻平穩條件：如果非週期馬爾科夫鏈的狀態轉移矩陣PP和概率分佈π(x)π(x)對於所有的i,ji,j滿足：

π(i)P(i,j)=π(j)P(j,i)π(i)P(i,j)=π(j)P(j,i)

　　　　則稱概率分佈π(x)π(x)是狀態轉移矩陣PP的平穩分佈。

　　　　在M-H取樣中我們通過引入接受率使細緻平穩條件滿足。現在我們換一個思路。

　　　　從二維的資料分佈開始，假設π(x1,x2)π(x1,x2)是一個二維聯合資料分佈，觀察第一個特徵維度相同的兩個點A(x(1)1,x(1)2)A(x1(1),x2(1))和B(x(1)1,x(2)2)B(x1(1),x2(2))，容易發現下面兩式成立：

π(x(1)1,x(1)2)π(x(2)2|x(1)1)=π(x(1)1)π(x(1)2|x(1)1)π(x(2)2|x(1)1)π(x1(1),x2(1))π(x2(2)|x1(1))=π(x1(1))π(x2(1)|x1(1))π(x2(2)|x1(1))

π(x(1)1,x(2)2)π(x(1)2|x(1)1)=π(x(1)1)π(x(2)2|x(1)1)π(x(1)2|x(1)1)π(x1(1),x2(2))π(x2(1)|x1(1))=π(x1(1))π(x2(2)|x1(1))π(x2(1)|x1(1))

　　　　由於兩式的右邊相等，因此我們有：

π(x(1)1,x(1)2)π(x(2)2|x(1)1)=π(x(1)1,x(2)2)π(x(1)2|x(1)1)π(x1(1),x2(1))π(x2(2)|x1(1))=π(x1(1),x2(2))π(x2(1)|x1(1))

　　　　也就是：

π(A)π(x(2)2|x(1)1)=π(B)π(x(1)2|x(1)1)π(A)π(x2(2)|x1(1))=π(B)π(x2(1)|x1(1))

　　　　觀察上式再觀察細緻平穩條件的公式，我們發現在x1=x(1)1x1=x1(1)這條直線上，如果用條件概率分佈π(x2|x(1)1)π(x2|x1(1))作為馬爾科夫鏈的狀態轉移概率，則任意兩個點之間的轉移滿足細緻平穩條件！這真是一個開心的發現，同樣的道理，在在x2=x(1)2x2=x2(1)這條直線上，如果用條件概率分佈π(x1|x(1)2)π(x1|x2(1))作為馬爾科夫鏈的狀態轉移概率，則任意兩個點之間的轉移也滿足細緻平穩條件。那是因為假如有一點C(x(2)1,x(1)2)C(x1(2),x2(1)),我們可以得到：

π(A)π(x(2)1|x(1)2)=π(C)π(x(1)1|x(1)2)π(A)π(x1(2)|x2(1))=π(C)π(x1(1)|x2(1))

　　　　基於上面的發現，我們可以這樣構造分佈π(x1,x2)π(x1,x2)的馬爾可夫鏈對應的狀態轉移矩陣PP：

P(A→B)=π(x(B)2|x(1)1)ifx(A)1=x(B)1=x(1)1P(A→B)=π(x2(B)|x1(1))ifx1(A)=x1(B)=x1(1)

P(A→C)=π(x(C)1|x(1)2)ifx(A)2=x(C)2=x(1)2P(A→C)=π(x1(C)|x2(1))ifx2(A)=x2(C)=x2(1)

P(A→D)=0elseP(A→D)=0else

　　　　有了上面這個狀態轉移矩陣，我們很容易驗證平面上的任意兩點E,FE,F，滿足細緻平穩條件：

π(E)P(E→F)=π(F)P(F→E)π(E)P(E→F)=π(F)P(F→E)

2.  二維Gibbs取樣

　　　　利用上一節找到的狀態轉移矩陣，我們就得到了二維Gibbs取樣，這個取樣需要兩個維度之間的條件概率。具體過程如下：

　　　　1）輸入平穩分佈π(x1,x2)π(x1,x2)，設定狀態轉移次數閾值n1n1，需要的樣本個數n2n2

　　　　2）隨機初始化初始狀態值x(0)1x1(0)和x(0)2x2(0)

　　　　3）for t=0t=0 to n1+n2−1n1+n2−1: 

　　　　　　a) 從條件概率分佈P(x2|x(t)1)P(x2|x1(t))中取樣得到樣本xt+12x2t+1

　　　　　　b) 從條件概率分佈P(x1|x(t+1)2)P(x1|x2(t+1))中取樣得到樣本xt+11x1t+1

　　　　樣本集{(x(n1)1,x(n1)2),(x(n1+1)1,x(n1+1)2),...,(x(n1+n2−1)1,x(n1+n2−1)2)}{(x1(n1),x2(n1)),(x1(n1+1),x2(n1+1)),...,(x1(n1+n2−1),x2(n1+n2−1))}即為我們需要的平穩分佈對應的樣本集。

　　　　整個取樣過程中，我們通過輪換座標軸，取樣的過程為：

(x(1)1,x(1)2)→(x(1)1,x(2)2)→(x(2)1,x(2)2)→...→(x(n1+n2−1)1,x(n1+n2−1)2)(x1(1),x2(1))→(x1(1),x2(2))→(x1(2),x2(2))→...→(x1(n1+n2−1),x2(n1+n2−1))

　　　　用下圖可以很直觀的看出，取樣是在兩個座標軸上不停的輪換的。當然，座標軸輪換不是必須的，我們也可以每次隨機選擇一個座標軸進行取樣。不過常用的Gibbs取樣的實現都是基於座標軸輪換的。

3. 多維Gibbs取樣

　　　　上面的這個演算法推廣到多維的時候也是成立的。比如一個n維的概率分佈π(x1,x2,...xn)π(x1,x2,...xn)，我們可以通過在n個座標軸上輪換取樣，來得到新的樣本。對於輪換到的任意一個座標軸xixi上的轉移，馬爾科夫鏈的狀態轉移概率為P(xi|x1,x2,...,xi−1,xi+1,...,xn)P(xi|x1,x2,...,xi−1,xi+1,...,xn)，即固定n−1n−1個座標軸，在某一個座標軸上移動。

　　　　具體的演算法過程如下：

　　　　1）輸入平穩分佈π(x1,x2，...,xn)π(x1,x2，...,xn)或者對應的所有特徵的條件概率分佈，設定狀態轉移次數閾值n1n1，需要的樣本個數n2n2

　　　　2）隨機初始化初始狀態值(x(0)1,x(0)2,...,x(0)n(x1(0),x2(0),...,xn(0)

　　　　3）for t=0t=0 to n1+n2−1n1+n2−1: 

　　　　　　a) 從條件概率分佈P(x1|x(t)2,x(t)3,...,x(t)n)P(x1|x2(t),x3(t),...,xn(t))中取樣得到樣本xt+11x1t+1

　　　　　　b) 從條件概率分佈P(x2|x(t+1)1,x(t)3,x(t)4,...,x(t)n)P(x2|x1(t+1),x3(t),x4(t),...,xn(t))中取樣得到樣本xt+12x2t+1

　　　　　　c)...

　　　　　　d) 從條件概率分佈P(xj|x(t+1)1,x(t+1)2,...,x(t+1)j−1,x(t)j+1...,x(t)n)P(xj|x1(t+1),x2(t+1),...,xj−1(t+1),xj+1(t)...,xn(t))中取樣得到樣本xt+1jxjt+1

　　　　　　e)...

　　　　　　f) 從條件概率分佈P(xn|x(t+1)1,x(t+1)2,...,x(t+1)n−1)P(xn|x1(t+1),x2(t+1),...,xn−1(t+1))中取樣得到樣本xt+1nxnt+1

　　　　樣本集{(x(n1)1,x(n1)2,...,x(n1)n),...,(x(n1+n2−1)1,x(n1+n2−1)2,...,x(n1+n2−1)n)}{(x1(n1),x2(n1),...,xn(n1)),...,(x1(n1+n2−1),x2(n1+n2−1),...,xn(n1+n2−1))}即為我們需要的平穩分佈對應的樣本集。

　　　　整個取樣過程和Lasso迴歸的座標軸下降法演算法非常類似，只不過Lasso迴歸是固定n−1n−1個特徵，對某一個特徵求極值。而Gibbs取樣是固定n−1n−1個特徵在某一個特徵取樣。

　　　　同樣的，輪換座標軸不是必須的，我們可以隨機選擇某一個座標軸進行狀態轉移，只不過常用的Gibbs取樣的實現都是基於座標軸輪換的。

4. 二維Gibbs取樣例項

　　　　這裡給出一個Gibbs取樣的例子。假設我們要取樣的是一個二維正態分佈Norm(μ,Σ)Norm(μ,Σ),其中：

μ=(μ1,μ2)=(5,−1)μ=(μ1,μ2)=(5,−1)

Σ=(σ21ρσ1σ2ρσ1σ2σ22)=(1114)Σ=(σ12ρσ1σ2ρσ1σ2σ22)=(1114)

　　　　而取樣過程中的需要的狀態轉移條件分佈為：

P(x1|x2)=Norm(μ1+ρσ1/σ2(x2−μ2),(1−ρ2)σ21)P(x1|x2)=Norm(μ1+ρσ1/σ2(x2−μ2),(1−ρ2)σ12)

P(x2|x1)=Norm(μ2+ρσ2/σ1(x1−μ1),(1−ρ2)σ22)P(x2|x1)=Norm(μ2+ρσ2/σ1(x1−μ1),(1−ρ2)σ22)

　　　　具體的程式碼如下：

from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal
samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])

def p_ygivenx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))

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系列就在這裡結束吧。