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