0.引子

書接前文，在采樣方法（一）中我們講到了拒絕采樣、重要性采樣一系列的蒙特卡洛采樣方法，但這些方法在高維空間時都會遇到一些問題，因為很難找到非常合適的可采樣Q分布，同時保證采樣效率以及精準(zhǔn)度。
本文將會介紹采樣方法中最重要的一族算法，MCMC（Markov Chain Monte Carlo），在之前我們的蒙特卡洛模擬都是按照如下公式進(jìn)行的：
$$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）馬爾可夫模型中我們就說到了馬爾可夫模型的馬爾可夫鏈，簡單來說就是滿足馬爾可夫假設(shè)
$$P(s_n|s_0,s_1,...,s_{n?1})=P(s_n|s_{n?1})$$
的狀態(tài)序列，由馬爾可夫過程（Markov Process）生成。
一個馬爾可夫過程由兩部分組成，一是狀態(tài)集合$\Omega$,二是轉(zhuǎn)移概率矩陣$T$。
其流程是：選擇一個初始的狀態(tài)分布${\pi }_0$，然后進(jìn)行狀態(tài)的轉(zhuǎn)移：
$${\pi }_n={\pi }_{n?1}T$$
得到的${\pi }_0,{\pi }_1,{\pi }_2...{\pi }_n$狀態(tài)分布序列。

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

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

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

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

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

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

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

2.MCMC

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

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

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

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

2.1.Metropolis

2.1.1原理描述

Metropolis算法算是MCMC的開山鼻祖了，它這里假設(shè)我們已經(jīng)有了一個狀態(tài)轉(zhuǎn)移概率函數(shù)T來表示轉(zhuǎn)移概率，T(a,b)表示從a轉(zhuǎn)移到b的概率(這里T的選擇我們稍后再說),顯然通常情況下一個T是不滿足細(xì)致平穩(wěn)條件的：
$$\pi (a)T(a,b)?=\pi (b)T(b,a)$$
所以我們需要進(jìn)行一些改造，加入一項$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}$$
所以我們改造后的滿足細(xì)致平穩(wěn)條件的轉(zhuǎn)移矩陣就是：
$$\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項是此次轉(zhuǎn)移的接受概率，是不是和拒絕采樣有點神似。

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

• 根據(jù)$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 })$，這個我們之前提到了是非常困難的，再回到我們的細(xì)致平穩(wěn)條件：
$$\pi (a)T(a,b)Q(a,b)=\pi (b)T(b,a)Q(b,a)$$
如果兩邊同時乘以一個數(shù)值，它也是成立的，比如
$$\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.有初始狀態(tài)$x_0$
2.從t=0,1,2,3…開始：

• 根據(jù)$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ài)轉(zhuǎn)移函數(shù)T的選擇，可以看到如果我們選擇一個對稱的轉(zhuǎn)移函數(shù)，即$T(a,b)=T(b,a)$，上面的接受概率還可以簡化為
$$\bar{Q}(a,b)=\min \left(1,\frac{\tilde{\pi }(b)}{\tilde{\pi }(a)}\right)$$
這也是一般Metropolis算法中采用的方法，T使用一個均勻分布即可，所有狀態(tài)之間的轉(zhuǎn)移概率都相同。
實驗中我們使用一個二元高斯分布來進(jìn)行采樣模擬

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

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

每輪采樣的函數(shù)：

def domain_random(): #計算定義域一個隨機(jī)值return np.random.random()*3.8-1.9 def metropolis(x):new_x = (domain_random(),domain_random()) #新狀態(tài)#計算接收概率acc = min(1,get_tilde_p((new_x[0],new_x[1]))/get_tilde_p((x[0],x[1])))#使用一個隨機(jī)數(shù)判斷是否接受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] #采樣狀態(tài)序列for i in range(counts):xs.append(x)x = metropolis(x) #采樣并判斷是否接受#在各個狀態(tài)之間繪制跳轉(zhuǎn)的線條幫助可視化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

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

2.2.Gibbs Sampling

2.2.1.算法分析

通過分析Metropolis的采樣軌跡，我們發(fā)現(xiàn)前后兩個狀態(tài)之間并沒有特別的聯(lián)系，新的狀態(tài)都是從T采樣出來的，而因為原始的分布很難計算，所以我們選擇的T是均勻分布，因此必須以一個概率進(jìn)行拒絕，才能保證最后收斂到我們期望的分布。
如果我們限定新的狀態(tài)只改變原狀態(tài)的其中一個維度，即
$$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})$為橋梁作轉(zhuǎn)換很好得到）：
$$\pi (x_t)?P(x_{new}^j|x_t^{?j})=\pi (x_{new})?P(x_t^j|x_t^{?j})$$

所以我們?nèi)绻麡?gòu)造這樣一個轉(zhuǎn)移概率函數(shù)T：
1.如果狀態(tài)a和b之間只在第j個元上值不同，則是的
$T(a,b)=P(b^{(j)}|a{(-j)}) $
2.否則$T(a,b)=0$

結(jié)論很清晰：這樣一個轉(zhuǎn)移概率函數(shù)T是滿足細(xì)致平穩(wěn)條件的，而且和Metropolis里面不同：它不是對稱的
我們能夠以1為概率接受它的轉(zhuǎn)移結(jié)果。

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

A只能跳轉(zhuǎn)到位于統(tǒng)一條坐標(biāo)線上的B,C兩個點，對于D，它無法一次轉(zhuǎn)移到達(dá)，但是可以通過兩次變換到達(dá)，仍然滿足Irreducible的條件。

這樣構(gòu)造出我們需要的轉(zhuǎn)移概率函數(shù)T之后，就直接得到我們的Gibbs采樣算法了：

Gibbs Sampling算法
1.有初始狀態(tài)$x_0$
2.從t=0,1,2,3…開始進(jìn)行循環(huán)采樣：

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

2.2.2.代碼實驗

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

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

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

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

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

def partialSampler(x,dim):xes = []for t in range(10): #隨機(jī)選擇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個點上進(jìn)行歸一化操作，然后按照概率進(jìn)行選擇。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]

主程序的結(jié)構(gòu)基本上和之前是一樣的，

def gibbs(x):rst = np.asarray(x)[:]path = [(x[0],x[1])]for dim in range(2): #維度輪詢，這里加入隨機(jī)也是可以的。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()

采樣的結(jié)果：

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

后注

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

• Metropolis是完全隨機(jī)的，以一個概率進(jìn)行拒絕；
• 而Gibbs Sampling則是在某個維度上進(jìn)行轉(zhuǎn)移。

如果我們?nèi)匀幌Ｍ詈笫褂?strong>獨立同分布的數(shù)據(jù)進(jìn)行蒙特卡洛模擬，只需要進(jìn)行多次MCMC，然后拿每個MCMC的第n個狀態(tài)作為一個樣本使用即可。

完整的代碼見鏈接：
https://github.com/justdark/dml/blob/master/dml/tool/mcmc.py
因為從頭到尾影響分布的只有get_p()函數(shù)，所以如果我們想對其他分布進(jìn)行實驗，只需要改變這個函數(shù)的定義就好了，比如說我們對一個相關(guān)系數(shù)為0.5的二元正態(tài)分布，只需要修改get_p()函數(shù)：
$$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))

就可以得到相應(yīng)的采樣結(jié)果：

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

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

也可以得到采樣結(jié)果：

PPS

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

Reference

【1】LDA數(shù)學(xué)八卦
【2】Pattern Recognition and Machine Learning
【3】Mathematicalmonk’s machine learning video

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。