0.引子

書接前文，在采樣方法（一）中我們講到了拒絕采樣、重要性采樣一系列的蒙特卡洛采樣方法，但這些方法在高維空間時都會遇到一些問題，因為很難找到非常合適的可采樣Q分布，同時保證采樣效率以及精準度。
本文將會介紹采樣方法中最重要的一族算法，MCMC（Markov Chain Monte Carlo），在之前我們的蒙特卡洛模擬都是按照如下公式進行的：
$$E[f(x)]\approx \frac{1}{m}\sum _{i=1}^{m}f(x_i).x_i～p.iid$$
我們的x都是獨立采樣出來的，而在MCMC中，它變成了
$$E[f(x)]\approx \frac{1}{m}\sum _{i=1}^{m}f(x_i).(x_0,x_1,...,x_m)～MC(p)$$
其中的$MC(p)$就是我們本文的主角之一，馬爾可夫過程（Markov Process）生成的馬爾可夫鏈（Markov Chain）。

1.Markov Chain（馬爾可夫鏈）

在序列的算法（一·a）馬爾可夫模型中我們就說到了馬爾可夫模型的馬爾可夫鏈，簡單來說就是滿足馬爾可夫假設
$$P(s_n|s_0,s_1,...,s_{n?1})=P(s_n|s_{n?1})$$
的狀態序列，由馬爾可夫過程（Markov Process）生成。
一個馬爾可夫過程由兩部分組成，一是狀態集合$\Omega$,二是轉移概率矩陣$T$。
其流程是：選擇一個初始的狀態分布${\pi }_0$，然后進行狀態的轉移：
$${\pi }_n={\pi }_{n?1}T$$
得到的${\pi }_0,{\pi }_1,{\pi }_2...{\pi }_n$狀態分布序列。

注意：我們在這里講的理論和推導都是基于離散變量的，但其實是可以直接推廣到連續變量。

馬爾可夫鏈在很多場景都有應用，比如大名鼎鼎的pagerank算法，都用到了類似的轉移過程；
馬爾可夫鏈在某種特定情況下，有一個奇妙的性質：

在某種條件下，當你從一個分布${\pi }_0$開始進行概率轉移，無論你一開始${\pi }_0$的選擇是什么，最終會收斂到一個固定的分布$\pi$，叫做穩態（steady-state）。
穩態滿足條件：
$$\pi =\pi T$$
這里可以參考《LDA數學八卦0.4.2》的例子，非常生動地描述了社會階層轉化的一個例子，也對MCMC作了非常好的講解

書歸正傳，回到我們采樣的場景，我們知道，采樣的難點就在于概率密度函數過于復雜而無法進行有效采樣，如果我們可以設計一個馬爾可夫過程，使得它最終收斂的分布是我們想要采樣的概率分布，不就可以解決我們的問題了么。

前面提到了在某種特定情況下，這就是所有MCMC算法的理論基礎Ergodic Theorem：
如果一個離散馬爾可夫鏈$(x_0,x_1...x_m)$是一個與時間無關的Irreducible的離散，并且有一個穩態分布$\pi$,則：
$$E[f(x)]\approx \frac{1}{m}\sum _{i=1}^{m}f(x_i).x～\pi$$
它需要滿足的條件有這樣幾個，我們直接列在這里，不作證明：

1.Time homogeneous: 狀態轉移與時間無關，這個很好理解。
2.Stationary Distribution: 最終是會收斂到穩定狀態的。
3.Irreducible: 任意兩個狀態之間都是可以互相到達的。
4.Aperiodic：馬爾可夫序列是非周期的，我們所見的絕大多數序列都是非周期的。

雖然這里要求是離散的馬爾可夫鏈，但實際上對于連續的場景也是適用的，只是轉移概率矩陣變成了轉移概率函數。

2.MCMC

在上面馬爾可夫鏈中我們的所說的狀態都是某個可選的變量值，比如社會等級上、中、下，而在采樣的場景中,特別是多元概率分布，并不是量從某個維度轉移到另一個維度，比如一個二元分布，二維平面上的每一個點都是一個狀態，所有狀態的概率和為1! 這里比較容易產生混淆，一定小心。

在這里我們再介紹一個概念：
Detail balance： 一個馬爾可夫過程是細致平穩的，即對任意a,b兩個狀態：
$$\pi (a)T_{ab}=\pi (b)T_{ba}$$細致平穩條件也可以推導出一個非周期的馬爾可夫鏈是平穩的，因為每次轉移狀態i從狀態j獲得的量與j從i獲得的量是一樣的，那毫無疑問最后$\pi T=\pi$.

所以我們的目標就是需要構造這樣一個馬爾可夫鏈，使得它最后能夠收斂到我們期望的分布$\pi$，而我們的狀態集合其實是固定的，所以最終目標就是構造一個合適的T，就大功告成了。

一般來說我們有:
$$\pi (x)=\frac{\tilde{\pi }(x)}{Z}$$
其中$Z$是歸一化參數$Z={\Sigma }_{x^{\prime }}\tilde{\pi }(x^{\prime })$，因為我們通常能夠很方便地計算分子$\tilde{\pi }(x)$，但是分母的計算因為要枚舉所有的狀態所以過于復雜而無法計算。我們希望最終采樣出來的樣本符合$\pi$分布。

2.1.Metropolis

2.1.1原理描述

Metropolis算法算是MCMC的開山鼻祖了，它這里假設我們已經有了一個狀態轉移概率函數T來表示轉移概率，T(a,b)表示從a轉移到b的概率(這里T的選擇我們稍后再說),顯然通常情況下一個T是不滿足細致平穩條件的：
$$\pi (a)T(a,b)?=\pi (b)T(b,a)$$
所以我們需要進行一些改造，加入一項$Q$使得等式成立：
$$\pi (a)T(a,b)Q(a,b)=\pi (b)T(b,a)Q(b,a)$$
基于對稱的原則，我們直接讓
$$\begin{gathered} Q(a,b)=\pi (b)T(b,a). \\ Q(b,a)=\pi (a)T(a,b). \end{gathered}$$
所以我們改造后的滿足細致平穩條件的轉移矩陣就是：
$$\begin{array}{rl}{T}^{\prime }(a,b)& =T(a,b)Q(a,b)\\ & =T(a,b)\left(\pi (b)T(b,a)\right)\end{array}$$
在Metropolis算法中,這個加入的這個Q項是此次轉移的接受概率，是不是和拒絕采樣有點神似。

MCMC采樣算法
1.有初始狀態$x_0$
2.從t=0,1,2,3開始：

• 根據$T(x|x_t)$采樣出$x^{\prime }$
• 以概率$Q(x_t,x^{\prime })=\pi (x^{\prime })T(x^{\prime },x_t)$接受，即$x_{t+1}=x^{\prime }$
• 否則使$x_{t+1}=x_t$

但這里還有一個問題，我們的接受概率Q可能會非常小，而且其中還需要精確計算出$\pi (x^{\prime })$，這個我們之前提到了是非常困難的，再回到我們的細致平穩條件：
$$\pi (a)T(a,b)Q(a,b)=\pi (b)T(b,a)Q(b,a)$$
如果兩邊同時乘以一個數值，它也是成立的，比如
$$\pi (a)T(a,b)Q(a,b)?10=\pi (b)T(b,a)Q(b,a)?10$$
所以我們可以同步放大Q(a,b)和Q(b,a)，使得其中最大的一個值為1，也就是說:
$$\begin{array}{rl}\bar{Q}(a,b)& =\min \left(1,\frac{Q(a,b)}{Q(b,a)}\right)\\ & =\min \left(1,\frac{\pi (b)T(b,a)}{\pi (a)T(a,b)}\right)\\ & =\min \left(1,\frac{\tilde{\pi }(b)T(b,a)}{\tilde{\pi }(a)T(a,b)}\right)\end{array}$$
這樣在提高接受率的同時，因為除式的存在我們還可以約掉難以計算的Z。

Metropolis-Hastings算法
1.有初始狀態$x_0$
2.從t=0,1,2,3…開始：

• 根據$T(x|x_t)$采樣出$x^{\prime }$
• 以概率$\bar{Q}(a,b)=\min \left(1,\frac{\tilde{\pi }(b)T(b,a)}{\tilde{\pi }(a)T(a,b)}\right)$接受，即$x_{t+1}=x^{\prime }$
• 否則使$x_{t+1}=x_t$

2.1.2.代碼實驗

我們之前提到狀態轉移函數T的選擇，可以看到如果我們選擇一個對稱的轉移函數，即$T(a,b)=T(b,a)$，上面的接受概率還可以簡化為
$$\bar{Q}(a,b)=\min \left(1,\frac{\tilde{\pi }(b)}{\tilde{\pi }(a)}\right)$$
這也是一般Metropolis算法中采用的方法，T使用一個均勻分布即可，所有狀態之間的轉移概率都相同。
實驗中我們使用一個二元高斯分布來進行采樣模擬

其概率密度函數這樣計算的,x是一個二維坐標：
$$p(x)=\frac{1}{2\pi }{e}^{?{x(0)}^2?{x(1)}^2}$$

def get_p(x):# 模擬pi函數return 1/(2*PI)*np.exp(- x[0]**2 - x[1]**2) def get_tilde_p(x):# 模擬不知道怎么計算Z的PI，20這個值對于外部采樣算法來說是未知的，對外只暴露這個函數結果return get_p(x)*20

每輪采樣的函數：

def domain_random(): #計算定義域一個隨機值return np.random.random()*3.8-1.9 def metropolis(x):new_x = (domain_random(),domain_random()) #新狀態#計算接收概率acc = min(1,get_tilde_p((new_x[0],new_x[1]))/get_tilde_p((x[0],x[1])))#使用一個隨機數判斷是否接受u = np.random.random()if u<acc:return new_xreturn x

然后就可以完整地跑一個實驗了：

def testMetropolis(counts = 100,drawPath = False):plotContour() #可視化#主要邏輯x = (domain_random(),domain_random()) #x0xs = [x] #采樣狀態序列for i in range(counts):xs.append(x)x = metropolis(x) #采樣并判斷是否接受#在各個狀態之間繪制跳轉的線條幫助可視化if drawPath: plt.plot(map(lambda x:x[0],xs),map(lambda x:x[1],xs),'k-',linewidth=0.5)##繪制采樣的點plt.scatter(map(lambda x:x[0],xs),map(lambda x:x[1],xs),c = 'g',marker='.')plt.show()pass

可以看到采樣結果并沒有想象的那么密集，因為雖然我們提高了接受率，但還是會拒絕掉很多點，所以即使采樣了5000次，繪制的點也沒有密布整個畫面。

2.2.Gibbs Sampling

2.2.1.算法分析

通過分析Metropolis的采樣軌跡，我們發現前后兩個狀態之間并沒有特別的聯系，新的狀態都是從T采樣出來的，而因為原始的分布很難計算，所以我們選擇的T是均勻分布，因此必須以一個概率進行拒絕，才能保證最后收斂到我們期望的分布。
如果我們限定新的狀態只改變原狀態的其中一個維度，即
$$x_{new}=(x_t^0,x_t^1,..x_{new}^j..,x_t^n)$$
只改變了其中第j個維度，則有$$\begin{gathered} \pi (x_t)=P(x_t^{?j})?P(x_t^j|x_t^{?j}) \\ \pi (x_{new})=P(x_t^{?j})?P(x_{new}^j|x_t^{?j}) \end{gathered}$$其中$x^{?j}$表示除了第j元其他的變量。
所以有（以$P(x_t^{?j})$為橋梁作轉換很好得到）：
$$\pi (x_t)?P(x_{new}^j|x_t^{?j})=\pi (x_{new})?P(x_t^j|x_t^{?j})$$

所以我們如果構造這樣一個轉移概率函數T：
1.如果狀態a和b之間只在第j個元上值不同，則是的
$T(a,b)=P(b^{(j)}|a{(-j)}) $
2.否則$T(a,b)=0$

結論很清晰：這樣一個轉移概率函數T是滿足細致平穩條件的，而且和Metropolis里面不同：它不是對稱的
我們能夠以1為概率接受它的轉移結果。

以一個二元分布為例，在平面上：

A只能跳轉到位于統一條坐標線上的B,C兩個點，對于D，它無法一次轉移到達，但是可以通過兩次變換到達，仍然滿足Irreducible的條件。

這樣構造出我們需要的轉移概率函數T之后，就直接得到我們的Gibbs采樣算法了：

Gibbs Sampling算法
1.有初始狀態$x_0$
2.從t=0,1,2,3…開始進行循環采樣：

• 將$x_t$復制過來，主要是為了循環采樣：$x_{t+1}=x_t$
• 枚舉維度j = 0…N，修改每個維度的值：
• • • 使得$x_{t+1}^j～P(x^j|x_{t+1}^{?j})$

2.2.2.代碼實驗

想必大家發現了，如果要寫代碼，我們必須要知道如何從$P(x^j|x_{t+1}^{?j})$進行采樣，這是一個后驗的概率分布，在很多應用中是能夠定義出函數表達的。

在我們之前實驗的場景（二元正態分布），確實也能精確地寫出這個概率分布的概率密度函數（也是一個正態分布）。

但退一步想，現在我們只關心一元的采樣了，所以其實是有很多方法可以用到的，比如拒絕采樣等等。

而最簡單的，直接在這一維度上隨機采幾個點，然后按照它們的概率密度函數值為權重選擇其中一個點作為采樣結果即可。

代碼里這樣做的目的主要是為了讓代碼足夠簡單，只依賴一個均勻分布的隨機數生成器。

def partialSampler(x,dim):xes = []for t in range(10): #隨機選擇10個點xes.append(domain_random())tilde_ps = []for t in range(10): #計算這10個點的未歸一化的概率密度值tmpx = x[:]tmpx[dim] = xes[t]tilde_ps.append(get_tilde_p(tmpx))#在這10個點上進行歸一化操作，然后按照概率進行選擇。norm_tilde_ps = np.asarray(tilde_ps)/sum(tilde_ps)u = np.random.random()sums = 0.0for t in range(10):sums += norm_tilde_ps[t]if sums>=u:return xes[t]

主程序的結構基本上和之前是一樣的，

def gibbs(x):rst = np.asarray(x)[:]path = [(x[0],x[1])]for dim in range(2): #維度輪詢，這里加入隨機也是可以的。new_value = partialSampler(rst,dim)rst[dim] = new_valuepath.append([rst[0],rst[1]])#這里最終只畫出一輪輪詢之后的點，但會把路徑都畫出來return rst,pathdef testGibbs(counts = 100,drawPath = False):plotContour()x = (domain_random(),domain_random())xs = [x]paths = [x]for i in range(counts):xs.append([x[0],x[1]])x,path = gibbs(x)paths.extend(path) #存儲路徑if drawPath:plt.plot(map(lambda x:x[0],paths),map(lambda x:x[1],paths),'k-',linewidth=0.5)plt.scatter(map(lambda x:x[0],xs),map(lambda x:x[1],xs),c = 'g',marker='.')plt.show()

采樣的結果：

其轉移的路徑看到都是與坐標軸平行的直線，并且可以看到采樣5000詞的時候跟Metropolis相比密集了很多，因為它沒有拒絕掉的點。

后注

本文我們講述了MCMC里面兩個最常見的算法Metropolis和Gibbs Sampling，以及它們各自的實現，從采樣路徑來看：

• Metropolis是完全隨機的，以一個概率進行拒絕；
• 而Gibbs Sampling則是在某個維度上進行轉移。

如果我們仍然希望最后使用獨立同分布的數據進行蒙特卡洛模擬，只需要進行多次MCMC，然后拿每個MCMC的第n個狀態作為一個樣本使用即可。

完整的代碼見鏈接：
https://github.com/justdark/dml/blob/master/dml/tool/mcmc.py
因為從頭到尾影響分布的只有get_p()函數，所以如果我們想對其他分布進行實驗，只需要改變這個函數的定義就好了，比如說我們對一個相關系數為0.5的二元正態分布，只需要修改get_p()函數：
$$p(x)=\frac{1}{2\pi \sqrt{1?0.5^2}}exp\left(?\frac{1}{2(1?0.5^2)}({x(0)}^2?x(0)x(1)+{x(1)}^2)\right)$$

def get_p(x):return 1/(2*PI*math.sqrt(1-0.25))*np.exp(-1/(2*(1-0.25))*(x[0]**2 -x[0]*x[1] + x[1]**2))

就可以得到相應的采樣結果：

而且因為這里并不要求p是一個歸一化后的分布，可以嘗試任何>0的函數，比如

def get_p(x):return np.sin(x[0]**2 + x[1]**2)+1

也可以得到采樣結果：

PPS

終于趕在2017年結束之前填了一個小坑，很久沒寫博客了。

Reference

【1】LDA數學八卦
【2】Pattern Recognition and Machine Learning
【3】Mathematicalmonk’s machine learning video

總結

以上是生活随笔為你收集整理的采样方法（二）MCMC相关算法介绍及代码实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。