MCMC隨機取樣
以下內容來自劉建平Pinard-部落格園的學習筆記，總結如下：
1 MCMC蒙特卡羅方法
作為一種隨機取樣方法，馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，以下簡稱MCMC）在機器學習，深度學習以及自然語言處理等領域都有廣泛的應用，是很多複雜演算法求解的基礎。下面我們就對MCMC的原理做一個總結。
1.1 MCMC概述
從名字我們可以看出，MCMC由兩個MC組成，即蒙特卡羅方法（Monte Carlo Simulation，簡稱MC）和馬爾科夫鏈（Markov Chain ，也簡稱MC）。要弄懂MCMC的原理我們首先得搞清楚蒙特卡羅方法和馬爾科夫鏈的原理。本節關注於蒙特卡羅方法。
1.2
蒙特卡羅方法引入
蒙特卡羅原來是一個賭場的名稱，用它作為名字大概是因為蒙特卡羅方法是一種隨機模擬的方法，這很像賭博場裡面的扔骰子的過程。最早的蒙特卡羅方法都是為了求解一些不太好求解的求和或者積分問題。比如積分：
如果我們很難求解出f（x）的原函式，那麼這個積分比較難求解。當然我們可以通過蒙特卡羅方法來模擬求解近似值。如何模擬呢？假設我們函式影象如下圖：
則一個簡單的近似求解方法是在［a，b］之間隨機的取樣一個點。比如
，然後用
代表在［a，b］區間上所有的
的值。那麼上面的定積分的近似求解為：
當然，用一個值代表［a，b］區間上所有的
的值，這個假設太粗糙。那麼我們可以取樣［a，b］區間的n個值：
，用它們的均值來代表［a，b］區間上所有的
的值。這樣我們上面的定積分的近似求解為：
雖然上面的方法可以一定程度上求解出近似的解，但是它隱含了一個假定，即x在［a，b］之間是均勻分佈的，而絕大部分情況，x在［a，b］之間不是均勻分佈的。如果我們用上面的方法，則模擬求出的結果很可能和真實值相差甚遠。
怎麼解決這個問題呢？
如果我們可以得到x在［a，b］的機率分佈函式p（x），那麼我們的定積分求和可以這樣進行：
上式最右邊的這個形式就是蒙特卡羅方法的一般形式。當然這裡是連續函式形式的蒙特卡羅方法，但是在離散時一樣成立。
可以看出，最上面我們假設x在［a，b］之間是均勻分佈的時候，
，帶入我們有機率分佈的蒙特卡羅積分的上式，可以得到：
也就是說，我們最上面的均勻分佈也可以作為一般機率分佈函式p（x）在均勻分佈時候的特例。那麼我們現在的問題轉到了如何求出x的分佈p（x）的若干和樣本上來。
1.3
機率分佈取樣
上一節我們講到蒙特卡羅方法的關鍵是得到x的機率分佈。如果求出了x的機率分佈，我們可以基於機率分佈去取樣基於這個機率分佈的n個x的樣本集，帶入蒙特卡羅求和的式子即可求解。但是還有一個關鍵的問題需要解決，即如何基於機率分佈去取樣基於這個機率分佈的n個x的樣本集。
其他一些常見的連續分佈，比如t分佈，F分佈，Beta分佈，Gamma分佈等，都可以透過類似的方式從uniform（0，1）得到的取樣樣本轉化得到。在python的numpy，scikit-learn等類庫中，都有生成這些常用分佈樣本的函式可以使用。
不過很多時候，我們的x的機率分佈不是常見的分佈，這意味著我們沒法方便的得到這些非常見的機率分佈的樣本集。那這個問題怎麼解決呢？
1.4
接受-拒絕取樣
對於機率分佈不是常見的分佈，一個可行的辦法是採用接受-拒絕取樣來得到該分佈的樣本。既然 p（x） 太複雜在程式中沒法直接取樣，那麼我設定一個程式可採樣的分佈 q（x） 比如高斯分佈，然後按照一定的方法拒絕某些樣本，以達到接近 p（x） 分佈的目的，其中q（x）叫做 proposal distribution。
具體採用過程如下，設定一個方便取樣的常用機率分佈函式 q（x），以及一個常量 k，使得 p（x） 總在 kq（x） 的下方。如上圖。
整個過程中，我們透過一系列的接受拒絕決策來達到用q（x）模擬p（x）機率分佈的目的。
1.5
蒙特卡羅方法小結
使用接受-拒絕取樣，我們可以解決一些機率分佈不是常見的分佈的時候，得到其取樣集並用蒙特卡羅方法求和的目的。但是接受-拒絕取樣也只能部分滿足我們的需求，在很多時候我們還是很難得到我們的機率分佈的樣本集。比如：
從上面可以看出，要想將蒙特卡羅方法作為一個通用的取樣模擬求和的方法，必須解決如何方便得到各種複雜機率分佈的對應的取樣樣本集的問題。馬爾科夫鏈就是幫助找到這些複雜機率分佈的對應的取樣樣本集的白衣騎士。
2 MCMC馬爾科夫鏈
用蒙特卡羅方法來隨機模擬求解一些複雜的連續積分或者離散求和的方法，但是這個方法需要得到對應的機率分佈的樣本集，而想得到這樣的樣本集很困難。因此我們需要馬爾科夫鏈來幫忙。
2.1 馬爾科夫鏈概述
馬爾科夫鏈定義本身比較簡單，它假設某一時刻狀態轉移的機率只依賴於它的前一個狀態。舉個形象的比喻，假如每天的天氣是一個狀態的話，那個今天是不是晴天只依賴於昨天的天氣，而和前天的天氣沒有任何關係。當然這麼說可能有些武斷，但是這樣做可以大大簡化模型的複雜度，因此馬爾科夫鏈在很多時間序列模型中得到廣泛的應用，比如迴圈神經網路RNN，隱式馬爾科夫模型HMM等，當然MCMC也需要它。
如果用精確的數學定義來描述，則假設我們的序列狀態是
既然某一時刻狀態轉移的機率只依賴於它的前一個狀態，那麼我們只要能求出系統中任意兩個狀態之間的轉換機率，這個馬爾科夫鏈的模型就定了。我們來看看下圖這個馬爾科夫鏈模型的具體的例子（來源於維基百科）。
這個馬爾科夫鏈是表示股市模型的，共有三種狀態：牛市（Bull market）， 熊市（Bear market）和橫盤（Stagnant market）。
每一個狀態都以一定的機率轉化到下一個狀態。比如，牛市以0。025的機率轉化到橫盤的狀態。這個狀態機率轉化圖可以以矩陣的形式表示。如果我們定義矩陣
某一位置
的值為
，即從狀態
轉化到狀態
的機率，並定義牛市為狀態0， 熊市為狀態1， 橫盤為狀態2。 這樣我們得到了馬爾科夫鏈模型的狀態轉移矩陣為：
馬爾科夫鏈模型的狀態轉移矩陣和我們蒙特卡羅方法需要的機率分佈樣本集有什麼關係呢？
這需要從馬爾科夫鏈模型的狀態轉移矩陣的性質講起。
2.2
馬爾科夫鏈模型狀態轉移矩陣的性質
得到了馬爾科夫鏈模型的狀態轉移矩陣，我們來看看馬爾科夫鏈模型的狀態轉移矩陣的性質。
仍然以上面的這個狀態轉移矩陣為例。假設我們當前股市的機率分佈為：［0。3，0。4，0。3］，即30%機率的牛市，40%機率的熊盤與30%的橫盤。然後這個狀態作為序列機率分佈的初始狀態
，將其帶入這個狀態轉移矩陣計算
的狀態。程式碼如下：
部分輸出結果如下：
Current round： 1
［［ 0。405 0。4175 0。1775］］
Current round： 2
［［ 0。4715 0。40875 0。11975］］
可以發現，從第60輪開始，我們的狀態機率分佈就不變了，一直保持在［0。625 0。3125 0。0625］，即62。5%的牛市，31。25%的熊市與6。25%的橫盤。那麼這個是巧合嗎？
我們現在換一個初始機率分佈試一試，現在我們用［0。7，0。1，0。2］作為初始機率分佈，然後這個狀態作為序列機率分佈的初始狀態
，將其帶入這個狀態轉移矩陣計算
的狀態。
程式碼如下：
部分輸出結果如下：
Current round： 1
［［ 0。695 0。1825 0。1225］］
Current round： 2
［［ 0。6835 0。22875 0。08775］］
可以看出，儘管這次我們採用了不同初始機率分佈，最終狀態的機率分佈趨於同一個穩定的機率分佈［0。625 0。3125 0。0625］， 也就是說我們的馬爾科夫鏈模型的狀態轉移矩陣收斂到的穩定機率分佈與我們的初始狀態機率分佈無關。這是一個非常好的性質，也就是說，如果我們得到了這個穩定機率分佈對應的馬爾科夫鏈模型的狀態轉移矩陣，則我們可以用任意的機率分佈樣本開始，帶入馬爾科夫鏈模型的狀態轉移矩陣，這樣經過一些序列的轉換，最終就可以得到符合對應穩定機率分佈的樣本。
這個性質不光對我們上面的狀態轉移矩陣有效，對於絕大多數的其他的馬爾科夫鏈模型的狀態轉移矩陣也有效。同時不光是離散狀態，連續狀態時也成立。
同時，對於一個確定的狀態轉移矩陣
，它的
次冪
在當n大於一定的值的時候也可以發現是確定的，我們還是以上面的例子為例，計算程式碼如下：
輸出結果如下：
Current round： 1
［［ 0。8275 0。13375 0。03875］
［ 0。2675 0。66375 0。06875］
［ 0。3875 0。34375 0。26875］］
Current round： 2
［［ 0。73555 0。212775 0。051675］
［ 0。42555 0。499975 0。074475］
［ 0。51675 0。372375 0。110875］］
我們可以發現，在n≥6以後，
的值穩定不再變化，而且每一行都為［0。625 0。3125 0。0625］，這和我們前面的穩定分佈是一致的。這個性質同樣不光是離散狀態，連續狀態時也成立。
用數學語言總結下馬爾科夫鏈的收斂性質：
上面的性質中需要解釋的有：
1）非週期的馬爾科夫鏈：這個主要是指馬爾科夫鏈的狀態轉化不是迴圈的，如果是迴圈的則永遠不會收斂。幸運的是我們遇到的馬爾科夫鏈一般都是非週期性的。用數學方式表述則是：對於任意某一狀態
，
為集合
的最大公約數，如果
，則該狀態為非週期的。
2）任何兩個狀態是連通的：這個指的是從任意一個狀態可以透過有限步到達其他的任意一個狀態，不會出現條件機率一直為0導致不可達的情況。
3）馬爾科夫鏈的狀態數可以是有限的，也可以是無限的。因此可以用於連續機率分佈和離散機率分佈。
4）
通常稱為馬爾科夫鏈的平穩分佈。
2.3
基於馬爾科夫鏈取樣
如果我們得到了某個平穩分佈所對應的馬爾科夫鏈狀態轉移矩陣，我們就很容易採用出這個平穩分佈的樣本集。
現在我們可以開始取樣了，首先，基於初始任意簡單機率分佈比如高斯分佈
取樣得到狀態值
，基於條件機率分佈
取樣狀態值
，一直進行下去，當狀態轉移進行到一定的次數時，比如到n次時，我們認為此時的取樣集
即是符合我們的平穩分佈的對應樣本集，可以用來做蒙特卡羅模擬求和了。
總結下基於馬爾科夫鏈的取樣過程：
2.4 馬爾科夫鏈取樣小結
如果假定我們可以得到我們需要取樣樣本的平穩分佈所對應的馬爾科夫鏈狀態轉移矩陣，那麼我們就可以用馬爾科夫鏈取樣得到我們需要的樣本集，進而進行蒙特卡羅模擬。
但是一個重要的問題是，隨意給定一個平穩分佈π，如何得到它所對應的馬爾科夫鏈狀態轉移矩陣P呢？
這是個大問題。
我們繞了一圈似乎還是沒有解決任意機率分佈取樣樣本集的問題。
幸運的是，MCMC取樣透過迂迴的方式解決了上面這個大問題，下面討論MCMC的取樣，以及它的使用改進版取樣： M-H取樣和Gibbs取樣。
3 MCMC取樣和M-H取樣
在MCMC馬爾科夫鏈中我們講到給定一個機率平穩分佈π， 很難直接找到對應的馬爾科夫鏈狀態轉移矩陣P。而只要解決這個問題，我們就可以找到一種通用的機率分佈取樣方法，進而用於蒙特卡羅模擬。本節我們就討論解決這個問題的辦法：MCMC取樣和它的易用版M-H取樣。
3.1 馬爾科夫鏈的細緻平穩條件
在解決從平穩分佈π， 找到對應的馬爾科夫鏈狀態轉移矩陣P之前，我們還需要先看看馬爾科夫鏈的細緻平穩條件。定義如下：
證明很簡單，由細緻平穩條件有：
將上式用矩陣表示即為：
即滿足馬爾可夫鏈的收斂性質。也就是說，只要我們找到了可以使機率分佈π（x）滿足細緻平穩分佈的矩陣P即可。這給了我們尋找從平穩分佈π， 找到對應的馬爾科夫鏈狀態轉移矩陣P的新思路。
不過不幸的是，僅僅從細緻平穩條件還是很難找到合適的矩陣P。比如我們的目標平穩分佈是π（x），隨機找一個馬爾科夫鏈狀態轉移矩陣Q，它是很難滿足細緻平穩條件的，即：
那麼如何使這個等式滿足呢？下面我們來看MCMC取樣如何解決這個問題。
3.2 MCMC取樣
由於一般情況下，目標平穩分佈π（x）和某一個馬爾科夫鏈狀態轉移矩陣Q不滿足細緻平穩條件，即
我們可以對上式做一個改造，使細緻平穩條件成立。方法是引入一個α（i，j），使上式可以取等號，即：
問題是什麼樣的α（i，j）可以使等式成立呢？其實很簡單，只要滿足下兩式即可：
這樣，就得到了分佈π（x）對應的馬爾科夫鏈狀態轉移矩陣P，滿足：
也就是說，我們的目標矩陣P可以透過任意一個馬爾科夫鏈狀態轉移矩陣Q乘以α（i，j）得到。
α（i，j）我們有一般稱之為接受率。取值在［0，1］之間，可以理解為一個機率值。即目標矩陣P可以透過任意一個馬爾科夫鏈狀態轉移矩陣Q以一定的接受率獲得。這個很像接受-拒絕取樣，那裡是以一個常用分佈透過一定的接受-拒絕機率得到一個非常見分佈，這裡是以一個常見的馬爾科夫鏈狀態轉移矩陣Q透過一定的接受-拒絕機率得到目標轉移矩陣P，兩者的解決問題思路是類似的。
總結下MCMC的取樣過程。
上面這個過程基本上就是MCMC取樣的完整取樣理論了，但是這個取樣演算法還是比較難在實際中應用，為什麼呢？問題在上面第三步的c步驟，接受率這兒。由於
可能非常的小，比如0。1，導致我們大部分的取樣值都被拒絕轉移，取樣效率很低。有可能我們取樣了上百萬次馬爾可夫鏈還沒有收斂，也就是上面這個
要非常非常的大，這讓人難以接受，怎麼辦呢？這時就輪到我們的M-H取樣出場了。
3.3
M-H取樣
M-H取樣是Metropolis-Hastings取樣的簡稱，這個演算法首先由Metropolis提出，被Hastings改進，因此被稱之為Metropolis-Hastings取樣或M-H取樣
M-H取樣解決了我們上一節MCMC取樣接受率過低的問題。
我們回到MCMC取樣的細緻平穩條件：
這時我們可以看到，如果兩邊同時擴大五倍，接受率提高到了0。5，但是細緻平穩條件卻仍然是滿足的，即：
這樣接受率可以做如下改進，即：
透過這個微小的改造，我們就得到了可以在實際應用中使用的M-H取樣演算法過程如下：
3.4
M-H取樣例項
為了更容易理解，這裡給出一個M-H取樣的例項。
在例子裡，我們的目標平穩分佈是一個均值3，標準差2的正態分佈，而選擇的馬爾可夫鏈狀態轉移矩陣Q（i，j）的條件轉移機率是以i為均值，方差1的正態分佈在位置j的值。這個例子僅僅用來讓大家加深對M-H取樣過程的理解。
程式碼如下：
輸出的圖中可以看到取樣值的分佈與真實的分佈之間的關係如下，取樣集還是比較擬合對應分佈的。
3.5
M-H取樣總結
M-H取樣完整解決了使用蒙特卡羅方法需要的任意機率分佈樣本集的問題，因此在實際生產環境得到了廣泛的應用。
但是在大資料時代，M-H取樣面臨著兩大難題：
2） 由於特徵維度大，很多時候我們甚至很難求出目標的各特徵維度聯合分佈，但是可以方便求出各個特徵之間的條件機率分佈。這時候我們能不能只有各維度之間條件機率分佈的情況下方便的取樣呢？
Gibbs取樣解決了上面兩個問題，因此在大資料時代，MCMC取樣基本是Gibbs取樣的天下。
4 MCMC:Gibbs取樣
M-H取樣已經可以很好的解決蒙特卡羅方法需要的任意機率分佈的樣本集的問題。
但是M-H取樣有兩個缺點：
一是需要計算接受率，在高維時計算量大。並且由於接受率的原因導致演算法收斂時間變長。
二是有些高維資料，特徵的條件機率分佈好求，但是特徵的聯合分佈不好求。
因此需要一個好的方法來改進M-H取樣，這就是我們下面講到的Gibbs取樣。
4.1
重新尋找合適的細緻平穩條件
細緻平穩條件：如果非週期馬爾科夫鏈的狀態轉移矩陣P和機率分佈π（x）對於所有的i，j滿足：
則稱機率分佈π（x）是狀態轉移矩陣P的平穩分佈。
在M-H取樣中我們透過引入接受率使細緻平穩條件滿足。現在我們換一個思路。
4.2
二維Gibbs取樣
利用上一節找到的狀態轉移矩陣，我們就得到了二維Gibbs取樣，這個取樣需要兩個維度之間的條件機率。具體過程如下：
用下圖可以很直觀的看出，取樣是在兩個座標軸上不停的輪換的。當然，座標軸輪換不是必須的，我們也可以每次隨機選擇一個座標軸進行取樣。不過常用的Gibbs取樣的實現都是基於座標軸輪換的。
4.3
多維Gibbs取樣
上面的這個演算法推廣到多維的時候也是成立的。比如一個n維的機率分佈
，可以透過在n個座標軸上輪換取樣，來得到新的樣本。對於輪換到的任意一個座標軸
上的轉移，馬爾科夫鏈的狀態轉移機率為
，即固定n−1個座標軸，在某一個座標軸上移動。
具體的演算法過程如下：
整個取樣過程和Lasso迴歸的座標軸下降法演算法非常類似，只不過Lasso迴歸是固定n−1個特徵，對某一個特徵求極值。而Gibbs取樣是固定n−1個特徵在某一個特徵取樣。
同樣的，輪換座標軸不是必須的，我們可以隨機選擇某一個座標軸進行狀態轉移，只不過常用的Gibbs取樣的實現都是基於座標軸輪換的。
4.4 二維Gibbs取樣例項
這裡給出一個Gibbs取樣的例子。假設我們要取樣的是一個二維正態分佈
，其中：
而取樣過程中的需要的狀態轉移條件分佈為：
具體的程式碼如下：
輸出的兩個特徵各自的分佈如下：
然後我們看看樣本集生成的二維正態分佈，程式碼如下：
輸出的正態分佈圖如下：
4.5
Gibbs取樣小結
由於Gibbs取樣在高維特徵時的優勢，目前我們通常意義上的MCMC取樣都是用的Gibbs取樣。
當然Gibbs取樣是從M-H取樣的基礎上的進化而來的，同時Gibbs取樣要求資料至少有兩個維度，一維機率分佈的取樣是沒法用Gibbs取樣的，這時M-H取樣仍然成立。
有了Gibbs取樣來獲取機率分佈的樣本集，有了蒙特卡羅方法來用樣本集模擬求和，他們一起就奠定了MCMC演算法在大資料時代高維資料模擬求和時的作用。