抽樣與蒙地卡羅(三)：馬可夫鏈蒙地卡羅方法

Sampling and Monte Carlo (03)： Markov Chain Monte Carlo (MCMC) Method

Markov Process and Markov Chain 馬可夫過程與馬可夫鏈介紹

在開始之前，先從基本的觀念複習起，記得在高中的時候有類似的經典問題(這邊參考周老師線代啟示錄部落格上的問題) :

如果台北每年有20%的人遷出外縣市，每年有10%的外縣市人口遷入台北市，假設總人口數一直不變，請問長期而言台北市與外縣市的人口比例會是如何 ? 針對這個問題，我們可以簡單列出一個公式， 令 s 代表外縣市人口數，t 代表台北市人口數，k 為年份，則可以有如下的方程解 :

或是改成矩陣表示 :

隨著 k 逼近無限大，我們可以解出最終該方程會趨於一個穩定比例，不論一開始外縣市與台北市人口比例是多少，最終都會趨於 2:1，此時我們將其稱為穩態(steady state)，詳細的計算請參考線代啟示錄，這邊就不多贅述。

馬可夫過程

當然，現實生活中要找出一個穩定的機率轉換矩陣幾乎是不太可能的，通常狀態轉換機率僅能通過近期狀態來推導出來，意思是也許2018年有20%的人遷出台北市，到了2019年，遷出的比例可能縮減成10%，也許是因為市內經濟發展提升，或是交通運輸方式改變等因素。因此，如果在狀態空間中，一個狀態轉移的隨機過程只與當前狀態有關，我們就稱之為馬可夫鏈(Markov Chain)，極端的例子就是隨機漫步，你永遠不可能預測下次狀態如何改變，僅能從當前的狀態去猜測可能的方向。

Source : Wikipedia ; 一個擁有兩個狀態的馬氏鏈

換算成數學形式則是 :

也就是給定過去 n 次狀態的情況下，第 n+1 次出現狀態的機率僅與上一次有關，至於其他性質我們下面再講，馬可夫鏈是個非常值得深入探討的題目。其應用也非常廣泛，舉凡本文要提到的積分與採樣方法，甚至到強化學習的發展，都與其有分不開的關係，而關於強化學習的一些雛形，我們留到下一篇再說，銜接上一篇，我們先回到積分問題，並引入馬可夫的概念。

馬可夫鏈性質與MCMC方法

上一篇的積分問題中，我們估計積分值的思路是從某特定分布中取樣，也就是將積分轉化成某機率密度下的期望值。在這篇我們則導入馬可夫鏈的概念，當馬可夫鏈收斂至平穩狀態 f(t) 時，可從其中抽樣來估計積分值。

整個問題在於， f(t) 有時候十分難以抽樣。此外，在大多數的情況下，觀察樣本與參數皆為隨機變量 :

這是比較完整的表達式。為甚麼要寫成這樣 ? 原因有幾個，其一、這樣的表達式能夠拓展到更高維度，很好地解決了我們上一篇難以事先從 f(t) 中取樣的問題。其二、這樣的表達式能夠計算邊際分布與後驗分布，對許多算法的推導上很有幫助。而這兩者都可以很好地透過一個叫做馬可夫鏈蒙地卡羅的方法來解決 (Markov Chain Monte Carlo，MCMC)，前者馬可夫鏈用以從目標分布中抽樣，後者蒙地卡羅用以模擬。

【馬可夫鏈性質與模擬 Properties and Simulation of Markov Chain】

而要回答原先的積分問題，我們首先知道馬可夫鏈會有最終將收斂到平穩分布的性質，因此只要確定了目標分布，我們可以將它作為平穩分布來創建馬可夫鏈，當運行足夠長的時間，該馬可夫鏈最終會收斂至該分布，也就可以用來抽樣並做積分近似了，而這就是整個 MCMC 方法的思路，但有個問題還沒解決，雖然馬可夫鏈的建構簡單，只要給定初始狀態與轉移矩陣，讓其不斷迭代就可以了，但怎麼能確保一開始的目標分布為我們給定的馬可夫鏈最終收斂的平穩分布，或反之如何根據已知的目標分布來設計模擬的馬可夫鏈 ? 這裡我們必須對馬可夫鏈做一些性質上的拆解 :

首先，馬可夫鏈如果反向來看，變成逆向的馬可夫鏈(reversed markov chain)，則可以用如下方式改寫 :

此時在時刻 n 的轉移機率變成了 :

以上，pi 定義為在時刻 n 時出現狀態 xn 的機率，會等於整條馬可夫鏈在時刻 n 時的邊際機率，而第一條公式是應用了貝氏定理中的條件機率性質來改寫，也就是反向來看從時刻 n+1 到時刻 n 的轉移機率，可以寫成時刻 n 到 n+1 的轉移機率乘上時刻 n 的邊際分布並除以時刻 n+1 的邊際分布，而因為我們知道當 n 足夠大時， pi 會收斂至平穩分布，因此有 :

這樣的方程式我們稱為 Detailed Balance，而得到這樣的性質以後，就很容易驗證目標分布是否為給定轉移矩陣下最終會收斂的平穩分布了，只要 :

其中最後一條等式即為平穩分布的定義，換個思路就是狀態等於 xn 的機率，等同於狀態 xn+1 的機率乘上從狀態 n+1 變成 n 的總機率 (逆向)。

有了這樣的式子，我們還是沒有辦法構建出馬可夫鏈，因為我們現在只是知道如何驗證我們提及的分佈在當前的馬可夫鏈下是最終的平穩分布，還是沒有提及我們如何「設計」一個符合收斂至該目標分布的馬可夫鏈。因此，這邊我們就來看一下，如何根據以上條件去設計演算法 :

【MCMC 算法步驟】

首先，如果隨機設計轉移矩陣 P，那麼很大機率不會符合上述的細緻平穩條件，也就是說 :

那麼，如果我們引入一個幫助用的算符，並假設等號成立，則 :

推導一下，令 :

其中 Q 我們稱為接受率 (Acceptance Rate)，數學意義上是用來幫助等式成立的。而直觀意義上，代表著當馬可夫鏈為平穩，則下一個狀態只與上一個狀態及其轉移機率有關；若非平穩，則下一個狀態與上一個狀態的轉移機率乘上接受率之後，可以視其為平穩序列並往下推進。因此，算法步驟為 :

1. 初始化起始狀態即轉移矩陣 P，其中轉移矩陣滿足對稱性，且總機率為 1；此外，定義須採樣的目標分布
2. 建構一條長度為 n 的馬可夫鏈 : 每次迭代遵照以下步驟 :
   一、針對當前狀態 i，從轉移矩陣 P 中得出提議狀態 j
   二、從均勻分布中取樣 u，如果 u < Q = pi * P，則轉移狀態至 j，否則停留在 i
   三、迭代 n 次，得出的馬可夫鏈即為目標分布中的採樣結果

可以來建立一個代碼，考慮一個簡單的離散目標分布，其狀態機率加總為 1，並隨機生成一個轉移矩陣，有四種不同狀態。此時我們來建構一個長度為 n的馬可夫鏈，也就是在該離散分布下採樣 n 個點，代碼如下 :

然而，此算法有一個缺點，當接受率以該方式定義時，容易造成過低而使拒絕率過高，此時容易因為初始狀態選擇或特定機率定義的問題而失準，讀者可以自己運行上面的程式碼並透過我們提過的驗證平穩分布的公式來嘗試檢驗看看，會發現十分不準，因此，以下提及一些改進的算法 :

Metropolis-Hastings 算法與抽樣

進入這裡之後，我們考慮狀態為連續的情況，前述的狀態都為離散的，計算相對簡單一些，但真實情況來說，很多時候要考慮連續狀態。

Metropolis-Hastings (以下簡稱 MH 算法)是上世紀 50 年代提出的經典算法，從基礎的 MCMC 上做出了一些改進，我們先看步驟 :

一、創造提議分布 g(x)；

二、從該提議分布 g(x) 裏頭產生候選狀態 J

三、若 J 被接受，則下一狀態為 J，否則停留在當前狀態 I

細心的讀者很快發現，這與剛才的 MCMC 方法不是完全一模一樣嗎 ? 沒錯，還記得我們剛才提到 MCMC 方法的缺陷在於容易出現接受率過低的問題，MH 算法實際上是針對這個接受率去做一些擴充。既然接受率過小容易造成停留在同一狀態，我們就對接收率做一個擴充 :

這樣的改法會限制接受率的最小為 1，避免了接受率無限下降的問題， min() 函數的右邊部分，考慮我們一開始談到的細緻平穩條件，我們可以知道右邊的等式最終為 1，而不斷輸出為 1 的接受率帶入我們上方說的擴充式，也可以反向確認滿足細緻平穩條件。

至此，我們某種程度上改善了接受率過低的問題，接下來另一個問題是，我們如何選擇好的提議分布 ? 雖然理論上來說，提議分布可以任意選取，並不會影響馬可夫鏈是否收斂至平穩分布，但考量計算效率，選擇一個與目標分布近似的分布效率會更好，通常以目標函數的支撐集(Support)為主，也就是集合 X 的一個子集，使在集合 X 上定義的目標函數 f 滿足在此子集上非 0。
** 至於為何選擇支撐集通常會有更好的效率，目前我不知道如何證明，參考了許多資料感覺是另一個數學領域，因此先留在這裡待補。

這裡我們透過程式碼實作練習看看，考慮一個目標分布為 :

我們選擇 sigma = 1，此外也選擇提議分布為 :

程式碼與模擬的結果如下 :

其中拒絕率約 6.9%。此外，讀著可以自行改變接受率的程式碼看看，會震盪到不要不要的XD

【Metropolis 方法中的 Weak Convergence and Optimal Scaling】

在上述的 MH 算法中，如果提議分布是對稱的，也就是 g(X|Y) = g(Y|X)的時候，我們稱此算法為 Metropolis 方法。

這邊要介紹的是在1997年時由 G.O. Roberts 等人提出，對普通 Metropolis 算法的一個補充，由於一般而言該算法十分仰賴自提議分布的伸縮，如果提議分布的方差過小，那麼收斂速度就會很慢，反之過大的話會造成拒絕率過高。一種改進的方法就是由外部監視接受率，當接受率在大約25%時，可以保證有比較好的收斂效果，因此該論文提議應該調整提議分布方法至運行時可以有近似的接受率。(文章提及的精確數字為 23%)

有興趣參考的讀者可以點進以下連結 :

Roberts , Gelman , Gilks : Weak convergence and optimal scaling of random walk Metropolis…

Gibbs 採樣算法 (Gibbs Sampler)

如果目標分布是多元分布，上述的MH算法就很難派上用場，因為推廣至多元分布後，計算上往往須考慮到難以計算的邊際分配；此外，提議分布的選取也會很大比例影響接受率。因此，Gibbs 採樣提出了一個基於各一元分布來採樣高維度的計算方式，而整個方法可以視作 MH 方法的一個特例。

核心概念如下，假設我們針對一個 n 維的分布做採樣，隨機變量為 X = (X1, X2, X3, …, Xn)，那麼我們可以定義任一個 n-1 維度的變量 :

吉布斯採樣的核心觀點在於從該條件分布中產生候選點，也就是改進了 MH 算法中一開始選擇的提議分布，讓他以後驗機率的形式與目標分布相連，這樣不僅改變了計算邊際分配的困難性，也讓提議分布與目標分布相關聯 :

一、針對 X 的各分量 1~n 迴圈，從提議分布 f(Xj | X-j) 中產生候選狀態 j，此外當前狀態為 i

二、更新狀態為 j (此時為針對所有分量迭代更新)

針對各分量做更新迭代的概念以圖呈現如下 :

Source : Little Orange’s Site

細心的讀者可以發現，這裡相比 MH 算法，大致步驟都一樣，唯獨缺少了一個接受率，在這裡我們可以將其看成接受率為 1，而至於為甚麼不需要擴充接受率的概念就可以確保滿足細緻平穩條件，實際上是因為將提議分布改成條件分布，在高維的情況下自然滿足此限制。

這裡的證明相對簡單，我就直接貼上卡內基美隆大學課程提供的等式了 :

Source : Carnegie Mellon University / Probabilistic Graphical Models

【Slice Gibbs Sampler 切片吉布斯抽樣】

當然，普通形式的吉布斯採樣並不完美，如果條件分布不好求，我們仍然難以透過程式實作，因此我們需要加入一些輔助用變量。

其中 u 必須要有方便抽樣的聯合分布，根據上面的算式，此時 x 的邊際分布就等於目標分布 f(x)，因此我們可以透過 g(u, x) 當作提議分布來產生隨機數，而當然選擇輔助變量非常重要，其聯合分布既然要好抽樣，一個簡單的選擇為均勻分布，此時 :

整個分布的撰寫就變得如此簡潔有力，一維圖解如下 :

這圖我拉好久…XD

總而言之，我們就可以嘗試來實作看看，假設我們目標產生一個二元的高斯分布 (python 有包可以調用，但這邊我們用 Gibbs 試試看):

Source : Wikipedia

當然，對於二元正態分布的參數，我們可以隨意定義，來比對一下我們透過切片抽樣與 Python Numpy 包中內建抽樣的結果。

左邊是透過 Slice Sampling 的結果，可以看到效果還是很不錯的。

結語

蒙地卡羅與抽樣方法到這裡告一段落。下一篇會進到這個系列的最後一章，但也可以說是番外篇，也就是關於隱馬可夫模型 (Hidden Markov Model，HMM) 的三大經典問題。該模型我認為是強化學習的入門磚，也是理解隨機過程中很重要的一環，雖然本質上並非取樣，但仍然沿襲我們這三篇所談的概念，因此我會將其放入番外篇介紹。

另外，最近幾次的文章都比較偏向數學面向，之後的隱馬可夫模型談完後，會更專注在應用面相關的文章，尤其是新穎的算法部分。