1. Introduction
第一次接觸到 Markov Chain Monte Carlo (MCMC) 是在 theano 的 deep learning tutorial 里面講解到的 RBM 用到了 Gibbs sampling，當(dāng)時(shí)因?yàn)橐s著做項(xiàng)目，雖然一頭霧水，但是也沒(méi)沒(méi)有時(shí)間仔細(xì)看。趁目前比較清閑，把 machine learning 里面的 sampling methods 理一理，發(fā)現(xiàn)內(nèi)容還真不少，有些知識(shí)本人也是一知半解，所以這篇AQF量化學(xué)習(xí)的文章不可能面面俱到詳細(xì)講解所有的 sampling methods，而是著重講一下這個(gè)號(hào)稱二十世紀(jì) top 10 之一的算法—— Markov chain Monte Carlo。在介紹 MCMC 之前，我們首先了解一下 MCMC 的 Motivation 和在它之前用到的方法。本人也是初學(xué)者，錯(cuò)誤在所難免，歡迎一起交流。
這篇文章從零開始，應(yīng)該都可以看懂，主要內(nèi)容包括：
隨機(jī)采樣
拒絕采樣
重要性采樣
Metropolis-Hastings Algorithm
Gibbs Sampling
2. Sampling
我們知道，計(jì)算機(jī)本身是無(wú)法產(chǎn)生真正的隨機(jī)數(shù)的，但是可以根據(jù)一定的算法產(chǎn)生偽隨機(jī)數(shù)(pseudo-random numbers)。最古老最簡(jiǎn)單的莫過(guò)于 Linear congruential generator：
式子中的 a 和 c 是一些數(shù)學(xué)知識(shí)推導(dǎo)出的合適的常數(shù)。但是我們看到，這種算法產(chǎn)生的下一個(gè)隨機(jī)數(shù)完全依賴現(xiàn)在的隨機(jī)數(shù)的大小，而且當(dāng)你的隨機(jī)數(shù)序列足夠大的時(shí)候，隨機(jī)數(shù)將出現(xiàn)重復(fù)子序列的情況。當(dāng)然，理論發(fā)展到今天，有很多更加先進(jìn)的隨機(jī)數(shù)產(chǎn)生算法出現(xiàn)，比如 python 數(shù)值運(yùn)算庫(kù) numpy 用的是 Mersenne Twister 等。但是不管算法如何發(fā)展，這些都不是本質(zhì)上的隨機(jī)數(shù)，用馮諾依曼的一句話說(shuō)就是：
Anyone who considers arithmetic methods of producing random digits is, of course, in a state of sin.
要檢查一個(gè)序列是否是真正的隨機(jī)序列，可以計(jì)算這個(gè)序列的 entropy 或者用壓縮算法計(jì)算該序列的冗余。
OK，根據(jù)上面的算法現(xiàn)在我們有了均勻分布的隨機(jī)數(shù)，但是如何產(chǎn)生滿足其他分布(比如高斯分布)下的隨機(jī)數(shù)呢?一種可選的簡(jiǎn)單的方法是 Inverse transform sampling，有時(shí)候也叫Smirnov transform。拿高斯分布舉例子，它的原理是利用高斯分布的累積分布函數(shù)(CDF，cumulative distribution function)來(lái)處理。
假如在 y 軸上產(chǎn)生(0,1)之間的均勻分布的隨機(jī)數(shù)，水平向右投影到高斯累計(jì)分布函數(shù)上，然后垂直向下投影到 x 軸，得到的就是高斯分布。可見高斯分布的隨機(jī)數(shù)實(shí)際就是均勻分布隨機(jī)數(shù)在高斯分布的 CDF 函數(shù)下的逆映射。當(dāng)然，在實(shí)際操作中，更有效的計(jì)算方法有 Box–Muller_transform (an efficient polar form)，Ziggurat algorithm 等，這些方法 tricky and faster，沒(méi)有深入了解，這里也不多說(shuō)了。
3. Motivation
MCMC 可解決高維空間里的積分和優(yōu)化問(wèn)題：
上面一個(gè)例子簡(jiǎn)單講了利用高斯分布的 CDF 可以產(chǎn)生高斯隨機(jī)數(shù)，但是有時(shí)候我們遇到一些分布的 CDF 計(jì)算不出來(lái)(無(wú)法用公式表示)，隨機(jī)數(shù)如何產(chǎn)生?
遇到某些無(wú)法直接求積分的函數(shù)，如 e^{x^2}，在計(jì)算機(jī)里面如何求積分?
如何對(duì)一個(gè)分布進(jìn)行高效快速的模擬，以便于抽樣?
如何在可行域很大(or large number of possible configurations)時(shí)有效找到較優(yōu)解——RBM 優(yōu)化目標(biāo)函數(shù)中的問(wèn)題。
如何在眾多模型中快速找到更好的模型——MDL, BIC, AIC 模型選擇問(wèn)題。
3.1 The Monte Carlo principle
實(shí)際上，Monte Carlo 抽樣基于這樣的思想：假設(shè)玩一局牌的贏的概率只取決于你抽到的牌，如果用窮舉的方法則有 52! 種情況，計(jì)算復(fù)雜度太大。而現(xiàn)實(shí)中的做法是先玩幾局試試，統(tǒng)計(jì)贏的概率，如果你不太確信這個(gè)概率，你可以盡可能多玩幾局，當(dāng)你玩的次數(shù)很大的時(shí)候，得到的概率就非常接近真實(shí)概率了。
上述方法可以估算隨機(jī)事件的概率，而用 Monte Carlo 抽樣計(jì)算隨即變量的期望值是接下來(lái)內(nèi)容的重點(diǎn)：X 表示隨即變量，服從概率分布 p(x), 那么要計(jì)算 f(x) 的期望，只需要我們不停從 p(x) 中抽樣
當(dāng)抽樣次數(shù)足夠的時(shí)候，就非常接近真實(shí)值了：
Monte Carlo 抽樣的方法還有一個(gè)重要的好處是：估計(jì)值的精度與 x 的維度無(wú)關(guān)(雖然維度越高，但是每次抽樣獲得的信息也越多)，而是與抽樣次數(shù)有關(guān)。在實(shí)際問(wèn)題里面抽樣二十次左右就能達(dá)到比較好的精度。
但是，當(dāng)我們實(shí)際動(dòng)手的時(shí)候，就會(huì)發(fā)現(xiàn)一個(gè)問(wèn)題——如何從分布 p(x) 中抽取隨機(jī)樣本。之前我們說(shuō)過(guò)，計(jì)算可以產(chǎn)生均勻分布的偽隨機(jī)數(shù)。顯然，第二小節(jié)產(chǎn)生高斯隨機(jī)數(shù)的抽樣方法只對(duì)少數(shù)特定的問(wèn)題管用，對(duì)于一般情況呢?
3.2 Rejection Sampling
Reject Sampling 實(shí)際采用的是一種迂回( proposal distribution q(x) )的策略。既然 p(x) 太復(fù)雜在程序中沒(méi)法直接采樣，那么我設(shè)定一個(gè)程序可抽樣的分布 q(x) 比如高斯分布，然后按照一定的方法拒絕某些樣本，達(dá)到接近 p(x) 分布的目的：
具體操作如下，設(shè)定一個(gè)方便抽樣的函數(shù) q(x)，以及一個(gè)常量 k，使得 p(x) 總在 kq(x) 的下方。(參考上圖)
x 軸方向：從 q(x) 分布抽樣得到 a。(如果是高斯，就用之前說(shuō)過(guò)的 tricky and faster 的算法更快)
y 軸方向：從均勻分布(0, kq(a)) 中抽樣得到 u。
如果剛好落到灰色區(qū)域： u > p(a), 拒絕， 否則接受這次抽樣
重復(fù)以上過(guò)程
在高維的情況下，Rejection Sampling 會(huì)出現(xiàn)兩個(gè)問(wèn)題，第一是合適的 q 分布比較難以找到，第二是很難確定一個(gè)合理的 k 值。這兩個(gè)問(wèn)題會(huì)導(dǎo)致拒絕率很高，無(wú)用計(jì)算增加。
3.3 Importance Sampling
Importance Sampling 也是借助了容易抽樣的分布 q (proposal distribution)來(lái)解決這個(gè)問(wèn)題，直接從公式出發(fā)：
其中，p(z) / q(z) 可以看做 importance weight。我們來(lái)考察一下上面的式子，p 和 f 是確定的，我們要確定的是 q。要確定一個(gè)什么樣的分布才會(huì)讓采樣的效果比較好呢?直觀的感覺(jué)是，樣本的方差越小期望收斂速率越快。比如一次采樣是 0, 一次采樣是 1000, 平均值是 500,這樣采樣效果很差，如果一次采樣是 499, 一次采樣是 501, 你說(shuō)期望是 500,可信度還比較高。在上式中，我們目標(biāo)是 p×f/q 方差越小越好，所以 |p×f| 大的地方，proposal distribution q(z) 也應(yīng)該大。舉個(gè)稍微極端的例子：
第一個(gè)圖表示 p 分布， 第二個(gè)圖的陰影區(qū)域 f = 1，非陰影區(qū)域 f = 0, 那么一個(gè)良好的 q 分布應(yīng)該在左邊箭頭所指的區(qū)域有很高的分布概率，因?yàn)樵谄渌麉^(qū)域的采樣計(jì)算實(shí)際上都是無(wú)效的。這表明 Importance Sampling 有可能比用原來(lái)的 p 分布抽樣更加有效。
但是可惜的是，在高維空間里找到一個(gè)這樣合適的 q 非常難。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出現(xiàn)，要找到一個(gè)同時(shí)滿足 easy to sample 并且 good approximations 的 proposal distribution, it is often impossible!
4. MCMC Algorithm
上面說(shuō)了這么多采樣方法，其實(shí)最終要突出的就是 MCMC 的過(guò)人之處。MCMC 的絕妙之處在于：通過(guò)穩(wěn)態(tài)的 Markov Chain 進(jìn)行轉(zhuǎn)移計(jì)算，等效于從 P(x) 分布采樣。但是在了解 MCMC 具體算法之前，我們還要熟悉 Markov Chain 是怎么一回事。
4.1 Markov Chain
Markov Chain 體現(xiàn)的是狀態(tài)空間的轉(zhuǎn)換關(guān)系，下一個(gè)狀態(tài)只決定與當(dāng)前的狀態(tài)(可以聯(lián)想網(wǎng)頁(yè)爬蟲原理，根據(jù)當(dāng)前頁(yè)面的超鏈接訪問(wèn)下一個(gè)網(wǎng)頁(yè))。如下圖：
這個(gè)狀態(tài)圖的轉(zhuǎn)換關(guān)系可以用一個(gè)轉(zhuǎn)換矩陣 T 來(lái)表示：
舉一個(gè)例子，如果當(dāng)前狀態(tài)為 u(x) = (0.5, 0.2, 0.3), 那么下一個(gè)矩陣的狀態(tài)就是 u(x)T = (0.18, 0.64, 0.18), 依照這個(gè)轉(zhuǎn)換矩陣一直轉(zhuǎn)換下去，最后的系統(tǒng)就趨近于一個(gè)穩(wěn)定狀態(tài) (0.22, 0.41, 0.37) (此處只保留了兩位有效數(shù)字)。而事實(shí)證明無(wú)論你從那個(gè)點(diǎn)出發(fā)，經(jīng)過(guò)很長(zhǎng)的 Markov Chain 之后都會(huì)匯集到這一點(diǎn)。
考慮一般的情況，滿足什么條件下經(jīng)過(guò)很長(zhǎng)的 Markov Chain 迭代后系統(tǒng)分布會(huì)趨近一個(gè)穩(wěn)定分布，即最后的 u(x) 等效于從目標(biāo)分布 p(x) 采樣。大概的條件如下(自己隨便總結(jié)的，可能有遺漏和錯(cuò)誤)：
Irreducibility. 即圖是聯(lián)通的，T 矩陣不能被切豆腐一樣劃分成小方塊，舉個(gè)例子，比如爬蟲爬不到內(nèi)部局域網(wǎng)的網(wǎng)頁(yè)
Aperiodicity. 即圖中遍歷不會(huì)陷入到一個(gè)死圈里，有些網(wǎng)站為了防機(jī)器人，會(huì)專門設(shè)置這種陷阱
Detailed Balance，這是保證系統(tǒng)有穩(wěn)態(tài)的一個(gè)重要條件，詳細(xì)說(shuō)明見下面。
假設(shè) p(x) 是最后的穩(wěn)態(tài)，那么 detailed balance 可以用公式表示為：
什么意思呢?假設(shè)上面狀態(tài)圖 x1 有 0.22 元， x2 有 0.41 元，x3 有 0.37 元，那么 0.22×1 表示 x1 需要給 x2 錢，以此類推，手動(dòng)計(jì)算，可以發(fā)現(xiàn)下一個(gè)狀態(tài)每個(gè)人手中的錢都沒(méi)有變。值得說(shuō)明的是，這里體現(xiàn)了一個(gè)很重要的特性，那就是從一個(gè)高概率狀態(tài) xi 向一個(gè)低概率狀態(tài) x(i-1) 轉(zhuǎn)移的概率等于從這個(gè)低概率狀態(tài)向高概率狀態(tài)轉(zhuǎn)移的概率(reversible，至于要不要轉(zhuǎn)移又是另外一回事)。當(dāng)然，在上面一個(gè)例子中，情況比較特殊，等號(hào)兩邊其實(shí)都是同一個(gè)東西。馬氏鏈的收斂性質(zhì)主要由轉(zhuǎn)移矩陣決定, 所以基于馬氏鏈做采樣的關(guān)鍵問(wèn)題是如何構(gòu)造轉(zhuǎn)移矩陣,使得平穩(wěn)分布恰好是我們要的分布p(x)。但是考慮一維的情況，假設(shè) p(x) 是一維高斯分布，x 是根據(jù) markov chain 得到的一個(gè)樣本，依照上面的等式，那么我們可以根據(jù)轉(zhuǎn)移矩陣 T左 和 T右 (這里實(shí)際是 proposal distribution)來(lái)得到 p(xi) 和 p(x(i-1)) 的比率，進(jìn)而按照一定的概率對(duì)這兩個(gè)樣本進(jìn)行選擇。通過(guò)大量這樣的處理，得到樣本就符合原始的 p(x) 分布了。這就是 MH 算法的基本原理。
4.2 Metropolis-Hastings Algorithm
舉個(gè)例子，我們要用 MH 算法對(duì)標(biāo)準(zhǔn)高斯分布進(jìn)行采樣，轉(zhuǎn)移函數(shù)(對(duì)稱)是方差為 0.05 的高斯矩陣，上述算法過(guò)程如下：
選取一個(gè)隨機(jī)點(diǎn) x0，作為一個(gè)采樣點(diǎn)
以 x0 為中心，以轉(zhuǎn)移函數(shù)為分布采取隨機(jī)點(diǎn) x1
以算法中的 A 概率接受 x1, 否則接受 x0
重復(fù)第二步第三步
注意到高斯分布是一個(gè)徑向基函數(shù)，上面算法畫波浪線的部分相等。
matlab 代碼如下：
n = 250000;x = zeros(n, 1);x(1) = 0.5;for i = 1: n-1 x_c = normrnd(x(1), 0.05); if rand < min(1, normpdf(x_c)/normpdf(x(i))) x(i+1) = x_c; else x(i+1) = x(i); end end
MH 算法中的 proposal distribution q(x) 也是需要小心確定的，詳細(xì)知識(shí)可以查閱這篇介紹論文 (An introduction to MCMC for machine learning, Andrieu, Christophe). 可以看到，這個(gè)算法和模擬退火算法的思想是非常相似的，但是在模擬退火算法過(guò)程中，隨著時(shí)間的增加，接受值大的區(qū)域的概率越來(lái)越高，直到找到較高點(diǎn)。
4.3 Gibbs Sampling
Gibbs Sampling 實(shí)際上是 MH 算法的一個(gè)變種。具體思路如下：假設(shè)在一定溫度下一定量的分子在容器里做無(wú)規(guī)則的熱運(yùn)動(dòng)，如何統(tǒng)計(jì)系統(tǒng)的能量呢?同樣，我們用 Monte Carlo 的思想進(jìn)行統(tǒng)計(jì)計(jì)算。我們假設(shè)所有的分子靜止在某一個(gè)時(shí)刻，這是初識(shí)狀態(tài)。固定其他的分子，根據(jù)分子間的作用力對(duì)其中一個(gè)分子進(jìn)行移動(dòng)，也就是說(shuō)在該分子以一定的概率移動(dòng)到領(lǐng)域的某一個(gè)地方，移動(dòng)完了之后再靜止。然后基于移動(dòng)后的狀態(tài)對(duì)下一個(gè)分子進(jìn)行同樣的移動(dòng)操作...直到所有的分子移動(dòng)完畢，那么現(xiàn)在的狀態(tài)就是 Monte Carlo 采樣的第二個(gè)樣本。依照這樣的順序采樣下去，我們對(duì)于這個(gè)系統(tǒng)就能計(jì)算一個(gè)統(tǒng)計(jì)意義上的能量了。從條件分布的角度來(lái)看，算法過(guò)程如下：
總體來(lái)講，Gibbs Sampling 就是從條件概率中選擇一個(gè)變量(分子)，然后對(duì)該變量(分子)進(jìn)行采樣。當(dāng)所有變量采樣完畢之后，就得到了后面的一個(gè)狀態(tài)，從而完成了對(duì)系統(tǒng)配置的采樣。在 deep learning 的 RBM 中，gibbs 采樣是已知權(quán)重參數(shù)和一個(gè) v 變量，通過(guò)采樣得到 h。通過(guò) h 采樣又可以得到另一個(gè) v ，如此交替采樣，就可以逐漸收斂于聯(lián)合分布了。
>>>點(diǎn)擊咨詢量化金融分析師AQF實(shí)訓(xùn)項(xiàng)目
.png)
完善下表，48小時(shí)內(nèi)查收全套AQF備考資料
.jpg)
