經(jīng)典的采樣方法有哪些？

▌引子

最近開始拾起來看一些 NLP 相關(guān)的東西，特別是深度學(xué)習(xí)在 NLP 上的應(yīng)用，發(fā)現(xiàn)采樣方法在很多模型中應(yīng)用得很多，因?yàn)橛?xùn)練的時(shí)候如果預(yù)測(cè)目標(biāo)是一個(gè)詞，直接的 softmax 計(jì)算量會(huì)根據(jù)單詞數(shù)量的增長而增長。恰好想到最開始深度學(xué)習(xí)在 DBN 的時(shí)候采樣也發(fā)揮了關(guān)鍵的作用，而自己對(duì)采樣相關(guān)的方法了解不算太多，所以去學(xué)習(xí)記錄一下，經(jīng)典的統(tǒng)計(jì)的方法確實(shí)巧妙，看起來非常有收獲。

本篇文章先主要介紹一下經(jīng)典的采樣方法如 Inverse Sampling、Rejective Sampling 以及 Importance Sampling 和它在 NLP 上的應(yīng)用，后面還會(huì)有一篇來嘗試介紹 MCMC 這一組狂炫酷拽的算法。才疏學(xué)淺，行文若有誤望指正。

▌Why Sampling

采樣是生活和機(jī)器學(xué)習(xí)算法中都會(huì)經(jīng)常用到的技術(shù)，一般來說采樣的目的是評(píng)估一個(gè)函數(shù)在某個(gè)分布上的期望值，也就是

比如我們都學(xué)過的拋硬幣，期望它的結(jié)果是符合一個(gè)伯努利分布的，定義正面的概率為 p ,反面概率為 1-p 。最簡單地使 f(x)=x，在現(xiàn)實(shí)中我們就會(huì)通過不斷地進(jìn)行拋硬幣這個(gè)動(dòng)作，來評(píng)估這個(gè)概率p。

這個(gè)方法也叫做蒙特卡洛法（Monte Carlo Method），常用于計(jì)算一些非常復(fù)雜無法直接求解的函數(shù)期望。

對(duì)于拋硬幣這個(gè)例子來說:

其期望就是拋到正面的計(jì)數(shù)除以總次數(shù) m。?

而我們拋硬幣的這個(gè)過程其實(shí)就是采樣，如果要用程序模擬上面這個(gè)過程也很簡單，因?yàn)椴植嫉臉颖竞苋菀咨桑?/p>

而在計(jì)算機(jī)中的隨機(jī)函數(shù)一般就是生成 0 到 1 的均勻分布隨機(jī)數(shù)。

▌Sampling Method

可以看到蒙特卡洛法其實(shí)就是按一定的概率分布中獲取大量樣本，用于計(jì)算函數(shù)在樣本的概率分布上的期望。其中最關(guān)鍵的一個(gè)步驟就是如何按照指定的概率分布 p 進(jìn)行樣本采樣，拋硬幣這個(gè) case 里伯努利分布是一個(gè)離散的概率分布，它的概率分布一般用概率質(zhì)量函數(shù)（pmf）表示，相對(duì)來說比較簡單，而對(duì)于連續(xù)概率分布我們需要考慮它的概率密度函數(shù)（pdf）：

比如上圖示例分別是標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)，它們的面積都是 1（這是概率的定義），如果我們可以按照相應(yīng)概率分布生成很多樣本，那這些樣本繪制出來的直方圖應(yīng)該跟概率密度函數(shù)是一致的。

而在實(shí)際的問題中，p 的概率密度函數(shù)可能會(huì)比較復(fù)雜，我們由淺入深，看看如何采樣方法如何獲得服從指定概率分布的樣本。

▌Inverse Sampling

對(duì)于一些特殊的概率分布函數(shù)，比如指數(shù)分布：

我們可以定義它的概率累積函數(shù)(Cumulative distribution function)，也就是(ps.這個(gè)’F’和前面的’f’函數(shù)并沒有關(guān)系)

從圖像上看就是概率密度函數(shù)小于 x 部分的面積。這個(gè)函數(shù)在 x≥0 的部分是一個(gè)單調(diào)遞增的函數(shù)(在定義域上單調(diào)非減)，定義域和值域是[0,+∞)→[0,1)，畫出來大概是這樣子的一個(gè)函數(shù)，在 p(x) 大的地方它增長快（梯度大），反之亦然：

因?yàn)樗俏ㄒ挥成涞模ㄔ?0的部分，接下來我們只考慮這一部分），所以它的反函數(shù)可以表示為，值域?yàn)閇0,+∞)

因?yàn)镕單調(diào)遞增，所以也是單調(diào)遞增的：?

利用反函數(shù)的定義，我們有：

我們定義一下 [0,1] 均勻分布的 CDF,這個(gè)很好理解：

所以

根據(jù) F(x) 的定義，它是 exp 分布的概率累積函數(shù)，所以上面這個(gè)公式的意思是符合 exp 分布，我們通過 F 的反函數(shù)將一個(gè) 0 到 1 均勻分布的隨機(jī)數(shù)轉(zhuǎn)換成了符合 exp 分布的隨機(jī)數(shù)，注意，以上推導(dǎo)對(duì)于 cdf 可逆的分布都是一樣的，對(duì)于 exp 來說，它的反函數(shù)的形式是：?

具體的映射關(guān)系可以看下圖(a)，我們從 y 軸 0-1 的均勻分布樣本（綠色）映射得到了服從指數(shù)分布的樣本（紅色）。

我們寫一點(diǎn)代碼來看看效果,最后繪制出來的直方圖可以看出來就是 exp 分布的圖，見上圖(b)，可以看到隨著采樣數(shù)量的變多，概率直方圖和真實(shí)的 CDF 就越接近：

defsampleExp(Lambda=2,maxCnt=50000):ys=[]standardXaxis=[]standardExp=[]foriinrange(maxCnt):u=np.random.random()y=-1/Lambda*np.log(1-u)#F-1(X)ys.append(y)foriinrange(1000):t=Lambda*np.exp(-Lambda*i/100)standardXaxis.append(i/100)standardExp.append(t)plt.plot(standardXaxis,standardExp,'r')plt.hist(ys,1000,normed=True)plt.show()

▌Rejective Sampling

我們?cè)趯W(xué)習(xí)隨機(jī)模擬的時(shí)候通常會(huì)講到用采樣的方法來計(jì)算 π 值，也就是在一個(gè) 1×1 的范圍內(nèi)隨機(jī)采樣一個(gè)點(diǎn)，如果它到原點(diǎn)的距離小于 1,則說明它在1/4 圓內(nèi)，則接受它，最后通過接受的占比來計(jì)算 1/4 圓形的面積，從而根據(jù)公式反算出預(yù)估的ππ值，隨著采樣點(diǎn)的增多，最后的結(jié)果會(huì)越精準(zhǔn)。?

上面這個(gè)例子里說明一個(gè)問題，我們想求一個(gè)空間里均勻分布的集合面積，可以嘗試在更大范圍內(nèi)按照均勻分布隨機(jī)采樣，如果采樣點(diǎn)在集合中，則接受，否則拒絕。最后的接受概率就是集合在‘更大范圍’的面積占比。

當(dāng)我們重新回過頭來看想要 sample 出來的樣本服從某一個(gè)分布 p，其實(shí)就是希望樣本在其概率密度函數(shù)高的地方出現(xiàn)得更多，所以一個(gè)直覺的想法，我們從均勻分布隨機(jī)生成一個(gè)樣本?，按照一個(gè)正比于?的概率接受這個(gè)樣本，也就是說雖然是從均勻分布隨機(jī)采樣，但留下的樣本更有可能是?高的樣本。

這樣的思路很自然，但是否是對(duì)的呢。其實(shí)這就是 Rejective Sampling 的基本思想，我們先看一個(gè)很 intuitive 的圖

假設(shè)目標(biāo)分布的 pdf 最高點(diǎn)是 1.5,有三個(gè)點(diǎn)它們的 pdf 值分別是

因?yàn)槲覀儚?x 軸上是按均勻分布隨機(jī)采樣的，所以采樣到三個(gè)點(diǎn)的概率都一樣，也就是

接下來需要決定每個(gè)點(diǎn)的接受概率，它應(yīng)該正比于?，當(dāng)然因?yàn)槭歉怕手狄残枰∮诘扔?1.?

我們可以畫一根 y=2 的直線，因?yàn)檎麄€(gè)概率密度函數(shù)都在這根直線下，我們?cè)O(shè)定

我們要做的就是生成一個(gè) 0-1的隨機(jī)數(shù)，如果它小于接受概率?，則留下這個(gè)樣本。因?yàn)?，所以可以看到因?yàn)?是?的3倍，所以?。同樣采集 100 次，最后留下來的樣本數(shù)期望也是 3 倍。這根本就是概率分布的定義!

我們將這個(gè)過程更加形式化一點(diǎn)，我們我們又需要采樣的概率密度函數(shù)

，但實(shí)際情況我們很有可能只能計(jì)算出?,有

我們需要找一個(gè)可以很方便進(jìn)行采樣的分布函數(shù)并使?

其中 c 是需要選擇的一個(gè)常數(shù)。然后我們從 q 分布中隨機(jī)采樣一個(gè)樣本,并以?

的概率決定是否接受這個(gè)樣本。重復(fù)這個(gè)過程就是「拒絕采樣」算法了。

在上面的例子我們選擇的 q 分布是均勻分布，所以從圖像上看其 pdf 是直線，但實(shí)際上和?越接近，采樣效率越高，因?yàn)槠浣邮芨怕室苍礁撸?

▌Importance Sampling

上面描述了兩種從另一個(gè)分布獲取指定分布的采樣樣本的算法，對(duì)于1.在實(shí)際工作中，一般來說我們需要 sample 的分布都及其復(fù)雜，不太可能求解出它的反函數(shù)，但 p(x) 的值也許還是可以計(jì)算的。對(duì)于2.找到一個(gè)合適的

那我們回過頭來看我們sample的目的：其實(shí)是想求得，也就是?

如果符合 p(x) 分布的樣本不太好生成，我們可以引入另一個(gè)分布 q(x)，可以很方便地生成樣本。使得

我們將問題轉(zhuǎn)化為了求 g(x) 在 q(x) 分布下的期望！?。?/p>

我們稱其中的叫做?Importance Weight.

▌Importance Sample 解決的問題

首先當(dāng)然是我們本來沒辦法 sample from p，這個(gè)是我們看到的，IS 將之轉(zhuǎn)化為了從 q 分布進(jìn)行采樣；同時(shí) IS 有時(shí)候還可以改進(jìn)原來的 sample，比如說:

可以看到如果我們直接從 p 進(jìn)行采樣，而實(shí)際上這些樣本對(duì)應(yīng)的 f(x) 都很小，采樣數(shù)量有限的情況下很有可能都無法獲得 f(x) 值較大的樣本，這樣評(píng)估出來的期望偏差會(huì)較大；

而如果我們找到一個(gè) q 分布，使得它能在 f(x)*p(x) 較大的地方采集到樣本，則能更好地逼近 [Ef(x)]，因?yàn)橛?Importance Weight 控制其比重，所以也不會(huì)導(dǎo)致結(jié)果出現(xiàn)過大偏差。

所以選擇一個(gè)好的p也能幫助你sample出來的效率更高，要使得 f(x)p(x)較大的地方能被 sample出來。

▌無法直接求得p(x)的情況

上面我們假設(shè) g(x) 和 q(x) 都可以比較方便地計(jì)算，但有些時(shí)候我們這個(gè)其實(shí)是很困難的，更常見的情況市我們能夠比較方便地計(jì)算和?

其中是一個(gè)標(biāo)準(zhǔn)化項(xiàng)（常數(shù)），使得?或者?等比例變化為一個(gè)概率分布，你可以理解為 softmax 里面那個(gè)除數(shù)。也就是說?

這種情況下我們的 importance sampling 是否還能應(yīng)用呢？

而我們直接計(jì)算并不太好計(jì)算，而它的倒數(shù)：?

因?yàn)槲覀兗以O(shè)能很方便地從 q 采樣，所以上式其實(shí)又被轉(zhuǎn)化成了一個(gè)蒙特卡洛可解的問題，也就是

最終最終，原來的蒙特卡洛問題變成了：

所以我們完全不用知道 q(x) 確切的計(jì)算值，就可以近似地從中得到在 q 分布下 f(x) 的取值??！amazing！

▌Importance Sampling在深度學(xué)習(xí)里面的應(yīng)用

在深度學(xué)習(xí)特別是NLP的Language Model中，訓(xùn)練的時(shí)候最后一層往往會(huì)使用 softmax 函數(shù)并計(jì)算相應(yīng)的梯度。

而我們知道 softmax 函數(shù)的表達(dá)式是：

要知道在 LM 中 m 的大小是詞匯的數(shù)量決定的，在一些巨大的模型里可能有幾十萬個(gè)詞，也就意味著計(jì)算Z的代價(jià)十分巨大。

而我們?cè)谟?xùn)練的時(shí)候無非是想對(duì) softmax 的結(jié)果進(jìn)行求導(dǎo)，也就是說

后面那一塊，我們好像看到了熟悉的東西，沒錯(cuò)這個(gè)形式就是為采樣量身定做似的。

經(jīng)典的蒙特卡洛方法就可以派上用途了，與其枚舉所有的詞，我們只需要從 V 里 sample 出一些樣本詞，就可以近似地逼近結(jié)果了。

同時(shí)直接從 P 中 sample 也不可取的，而且計(jì)算 P是非常耗時(shí)的事情（因?yàn)樾枰?jì)算Z），我們一般只能計(jì)算，而且直接從 P 中 sample 也不可取，所以我們選擇另一個(gè)分布 Q 進(jìn)行 Importance Sample 即可。

一般來說可能選擇的 Q 分布是簡單一些的 n-gramn 模型。下面是論文中的算法偽代碼，基本上是比較標(biāo)準(zhǔn)的流程（論文圖片的符號(hào)和上面的描述稍有出入，理解一下過程即可）：

References

【1】mathematicalmonk’s machine learning course on y2b. machine learing

【2】Pattern Recognition And Machine Learning

【3】Adaptive Importance Sampling to Accelerate Trainingof a Neural Probabilistic Language Model.Yoshua Bengio and Jean-Sébastien Senécal.

?

采樣方法 2

▌引子

在上面我們講到了拒絕采樣、重要性采樣一系列的蒙特卡洛采樣方法，但這些方法在高維空間時(shí)都會(huì)遇到一些問題，因?yàn)楹茈y找到非常合適的可采樣Q分布，同時(shí)保證采樣效率以及精準(zhǔn)度。

本文將會(huì)介紹采樣方法中最重要的一族算法，MCMC（Markov Chain Monte Carlo），在之前我們的蒙特卡洛模擬都是按照如下公式進(jìn)行的：

我們的x都是獨(dú)立采樣出來的，而在MCMC中，它變成了

其中的MC(p)就是我們本文的主角之一，馬爾可夫過程（Markov Process）生成的馬爾可夫鏈（Markov Chain）。

▌Markov Chain（馬爾可夫鏈）

在序列的算法（一·a）馬爾可夫模型中（https://blog.csdn.net/dark_scope/article/details/61417336）我們就說到了馬爾可夫模型的馬爾可夫鏈，簡單來說就是滿足馬爾可夫假設(shè)

的狀態(tài)序列，由馬爾可夫過程（Markov Process）生成。

一個(gè)馬爾可夫過程由兩部分組成，一是狀態(tài)集合Ω，二是轉(zhuǎn)移概率矩陣T。

其流程是：選擇一個(gè)初始的狀態(tài)分布π0，然后進(jìn)行狀態(tài)的轉(zhuǎn)移：

得到的π0，π1，π2 ...πn 狀態(tài)分布序列。

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

馬爾可夫鏈在很多場(chǎng)景都有應(yīng)用，比如大名鼎鼎的 pagerank 算法，都用到了類似的轉(zhuǎn)移過程；

馬爾可夫鏈在某種特定情況下，有一個(gè)奇妙的性質(zhì)：在某種條件下，當(dāng)你從一個(gè)分布π0開始進(jìn)行概率轉(zhuǎn)移，無論你一開始π0 的選擇是什么，最終會(huì)收斂到一個(gè)固定的分布π，叫做穩(wěn)態(tài)（steady-state）。

穩(wěn)態(tài)滿足條件：

這里可以參考《LDA數(shù)學(xué)八卦0.4.2》的例子，非常生動(dòng)地描述了社會(huì)階層轉(zhuǎn)化的一個(gè)例子，也對(duì)MCMC作了非常好的講解

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

前面提到了在某種特定情況下，這就是所有MCMC算法的理論基礎(chǔ)ErgodicTheorem：

如果一個(gè)離散馬爾可夫鏈(x0,x1...xm)是一個(gè)與時(shí)間無關(guān)的 Irreducible 的離散，并且有一個(gè)穩(wěn)態(tài)分布π,則：

它需要滿足的條件有這樣幾個(gè)，我們直接列在這里，不作證明：????????????????

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

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

▌MCMC

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

在這里我們?cè)俳榻B一個(gè)概念：

Detail balance：一個(gè)馬爾可夫過程是細(xì)致平穩(wěn)的，即對(duì)任意 a,b 兩個(gè)狀態(tài)：

細(xì)致平穩(wěn)條件也可以推導(dǎo)出一個(gè)非周期的馬爾可夫鏈?zhǔn)瞧椒€(wěn)的，因?yàn)槊看无D(zhuǎn)移狀態(tài)i從狀態(tài)j獲得的量與 j 從 i 獲得的量是一樣的，那毫無疑問最后πT=π.

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

一般來說我們有:

其中Z是歸一化參數(shù),因?yàn)槲覀兺ǔＤ軌蚝芊奖愕赜?jì)算分子,但是分母的計(jì)算因?yàn)橐杜e所有的狀態(tài)所以過于復(fù)雜而無法計(jì)算。我們希望最終采樣出來的樣本符合?π 分布。

▌Metropolis

原理描述

Metropolis 算法算是 MCMC 的開山鼻祖了，它這里假設(shè)我們已經(jīng)有了一個(gè)狀態(tài)轉(zhuǎn)移概率函數(shù)T來表示轉(zhuǎn)移概率，T(a,b) 表示從 a 轉(zhuǎn)移到 b 的概率(這里T的選擇我們稍后再說),顯然通常情況下一個(gè)T是不滿足細(xì)致平穩(wěn)條件的：

所以我們需要進(jìn)行一些改造，加入一項(xiàng)Q使得等式成立：

基于對(duì)稱的原則，我們直接讓

所以我們改造后的滿足細(xì)致平穩(wěn)條件的轉(zhuǎn)移矩陣就是：

在 Metropolis 算法中,這個(gè)加入的這個(gè) Q 項(xiàng)是此次轉(zhuǎn)移的接受概率，是不是和拒絕采樣有點(diǎn)神似。

但這里還有一個(gè)問題，我們的接受概率 Q 可能會(huì)非常小，而且其中還需要精確計(jì)算出π(x′)，這個(gè)我們之前提到了是非常困難的，再回到我們的細(xì)致平穩(wěn)條件：

如果兩邊同時(shí)乘以一個(gè)數(shù)值，它也是成立的，比如

所以我們可以同步放大Q(a,b)和Q(b,a)，使得其中最大的一個(gè)值為1，也就是說:

這樣在提高接受率的同時(shí)，因?yàn)槌降拇嬖谖覀冞€可以約掉難以計(jì)算的 Z。

代碼實(shí)驗(yàn)

我們之前提到狀態(tài)轉(zhuǎn)移函數(shù)T的選擇，可以看到如果我們選擇一個(gè)對(duì)稱的轉(zhuǎn)移函數(shù),即T(a,b)=T(b,a),上面的接受概率還可以簡化為

這也是一般 Metropolis 算法中采用的方法，T使用一個(gè)均勻分布即可，所有狀態(tài)之間的轉(zhuǎn)移概率都相同。

實(shí)驗(yàn)中我們使用一個(gè)二元高斯分布來進(jìn)行采樣模擬

其概率密度函數(shù)這樣計(jì)算的,x是一個(gè)二維坐標(biāo)：

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

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

defdomain_random():#計(jì)算定義域一個(gè)隨機(jī)值returnnp.random.random()*3.8-1.9defmetropolis(x):new_x=(domain_random(),domain_random())#新狀態(tài)#計(jì)算接收概率acc=min(1,get_tilde_p((new_x[0],new_x[1]))/get_tilde_p((x[0],x[1])))#使用一個(gè)隨機(jī)數(shù)判斷是否接受u=np.random.random()ifu

然后就可以完整地跑一個(gè)實(shí)驗(yàn)了：

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

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

▌Gibbs Sampling

算法分析

通過分析 Metropolis 的采樣軌跡，我們發(fā)現(xiàn)前后兩個(gè)狀態(tài)之間并沒有特別的聯(lián)系，新的狀態(tài)都是從 T 采樣出來的，而因?yàn)樵嫉姆植己茈y計(jì)算，所以我們選擇的T是均勻分布，因此必須以一個(gè)概率進(jìn)行拒絕，才能保證最后收斂到我們期望的分布。

如果我們限定新的狀態(tài)只改變?cè)瓲顟B(tài)的其中一個(gè)維度，即

只改變了其中第 j 個(gè)維度，則有

其中表示除了第j元其他的變量。?所以有（以為橋梁作轉(zhuǎn)換很好得到）：?

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

我們能夠以 1 為概率接受它的轉(zhuǎn)移結(jié)果。

以一個(gè)二元分布為例，在平面上：

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

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

代碼實(shí)驗(yàn)

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

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

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

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

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

defpartialSampler(x,dim):xes=[]fortinrange(10):#隨機(jī)選擇10個(gè)點(diǎn)xes.append(domain_random())tilde_ps=[]fortinrange(10):#計(jì)算這10個(gè)點(diǎn)的未歸一化的概率密度值tmpx=x[:]tmpx[dim]=xes[t]tilde_ps.append(get_tilde_p(tmpx))#在這10個(gè)點(diǎn)上進(jìn)行歸一化操作，然后按照概率進(jìn)行選擇。norm_tilde_ps=np.asarray(tilde_ps)/sum(tilde_ps)u=np.random.random()sums=0.0fortinrange(10):sums+=norm_tilde_ps[t]ifsums>=u:returnxes[t]

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

defgibbs(x):rst=np.asarray(x)[:]path=[(x[0],x[1])]fordiminrange(2):#維度輪詢，這里加入隨機(jī)也是可以的。new_value=partialSampler(rst,dim)rst[dim]=new_valuepath.append([rst[0],rst[1]])#這里最終只畫出一輪輪詢之后的點(diǎn)，但會(huì)把路徑都畫出來returnrst,pathdeftestGibbs(counts=100,drawPath=False):plotContour()x=(domain_random(),domain_random())xs=[x]paths=[x]foriinrange(counts):xs.append([x[0],x[1]])x,path=gibbs(x)paths.extend(path)#存儲(chǔ)路徑ifdrawPath:plt.plot(map(lambdax:x[0],paths),map(lambdax:x[1],paths),'k-',linewidth=0.5)plt.scatter(map(lambdax:x[0],xs),map(lambdax:x[1],xs),c='g',marker='.')plt.show()

采樣的結(jié)果：

??????????????

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

▌后注

本文我們講述了MCMC里面兩個(gè)最常見的算法Metropolis和Gibbs Sampling，以及它們各自的實(shí)現(xiàn)

從采樣路徑來看：

Metropolis 是完全隨機(jī)的，以一個(gè)概率進(jìn)行拒絕；

而 Gibbs Sampling 則是在某個(gè)維度上進(jìn)行轉(zhuǎn)移。

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

完整的代碼見鏈接：

https://github.com/justdark/dml/blob/master/dml/tool/mcmc.py

因?yàn)閺念^到尾影響分布的只有g(shù)et_p()函數(shù)，所以如果我們想對(duì)其他分布進(jìn)行實(shí)驗(yàn)，只需要改變這個(gè)函數(shù)的定義就好了，比如說我們對(duì)一個(gè)相關(guān)系數(shù)為0.5 的二元正態(tài)分布，只需要修改 get_p() 函數(shù)：

defget_p(x):return1/(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é)果：

而且因?yàn)檫@里并不要求p是一個(gè)歸一化后的分布，可以嘗試任何>0的函數(shù)，比如

defget_p(x):returnnp.sin(x[0]**2+x[1]**2)+1

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