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