MCMC(二)馬爾科夫鏈
在MCMC(一)蒙特卡羅方法中，我們講到了如何用蒙特卡羅方法來隨機模擬求解一些複雜的連續積分或者離散求和的方法，但是這個方法需要得到對應的概率分佈的樣本集，而想得到這樣的樣本集很困難。因此我們需要本篇講到的馬爾科夫鏈來幫忙。
1. 馬爾科夫鏈概述
馬爾科夫鏈定義本身比較簡單，它假設某一時刻狀態轉移的概率只依賴於它的前一個狀態。舉個形象的比喻，假如每天的天氣是一個狀態的話，那個今天是不是晴天只依賴於昨天的天氣，而和前天的天氣沒有任何關係。當然這麼說可能有些武斷，但是這樣做可以大大簡化模型的複雜度，因此馬爾科夫鏈在很多時間序列模型中得到廣泛的應用，比如迴圈神經網路RNN，隱式馬爾科夫模型HMM等，當然MCMC也需要它。
如果用精確的數學定義來描述，則假設我們的序列狀態是$...X_{t-2}, X_{t-1}, X_{t}, X_{t+1},...$，那麼我們的在時刻$X_{t+1}$的狀態的條件概率僅僅依賴於時刻$X_{t}$，即：$$P(X_{t+1} |...X_{t-2}, X_{t-1}, X_{t} ) = P(X_{t+1} | X_{t})$$
既然某一時刻狀態轉移的概率只依賴於它的前一個狀態，那麼我們只要能求出系統中任意兩個狀態之間的轉換概率，這個馬爾科夫鏈的模型就定了。我們來看看下圖這個馬爾科夫鏈模型的具體的例子(來源於維基百科)。
這個馬爾科夫鏈是表示股市模型的，共有三種狀態：牛市（Bull market）, 熊市（Bear market）和橫盤（Stagnant market）。每一個狀態都以一定的概率轉化到下一個狀態。比如，牛市以0.025的概率轉化到橫盤的狀態。這個狀態概率轉化圖可以以矩陣的形式表示。如果我們定義矩陣陣$P$某一位置$P(i,j)$的值為$P(j|i)$,即從狀態i轉化到狀態j的概率，並定義牛市為狀態0， 熊市為狀態1, 橫盤為狀態2. 這樣我們得到了馬爾科夫鏈模型的狀態轉移矩陣為：$$P=\left( \begin{array}{ccc}
0.9&0.075&0.025 \\ 0.15&0.8& 0.05 \\
0.25&0.25&0.5 \end{array} \right)$$
講了這麼多，那麼馬爾科夫鏈模型的狀態轉移矩陣和我們蒙特卡羅方法需要的概率分佈樣本集有什麼關係呢？這需要從馬爾科夫鏈模型的狀態轉移矩陣的性質講起。
2. 馬爾科夫鏈模型狀態轉移矩陣的性質
得到了馬爾科夫鏈模型的狀態轉移矩陣，我們來看看馬爾科夫鏈模型的狀態轉移矩陣的性質。
完整程式碼參見我的github:https://github.com/ljpzzz/machinelearning/blob/master/mathematics/mcmc_2.ipynb
仍然以上面的這個狀態轉移矩陣為例。假設我們當前股市的概率分佈為：[0.3,0.4,0.3],即30%概率的牛市，40%概率的熊盤與30%的橫盤。然後這個狀態作為序列概率分佈的初始狀態$t_0$，將其帶入這個狀態轉移矩陣計算$t_1, t_2,t_3...$的狀態。程式碼如下：
import numpy as np matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) vector1 = np.matrix([[0.3,0.4,0.3]], dtype=float) for i in range(100): vector1 = vector1*matrix print "Current round:" , i+1 print vector1
部分輸出結果如下：
Current round: 1 [[ 0.405 0.4175 0.1775]] Current round: 2 [[ 0.4715 0.40875 0.11975]] Current round: 3 [[ 0.5156 0.3923 0.0921]] Current round: 4 [[ 0.54591 0.375535 0.078555]] 。。。。。。 Current round: 58 [[ 0.62499999 0.31250001 0.0625 ]] Current round: 59 [[ 0.62499999 0.3125 0.0625 ]] Current round: 60 [[ 0.625 0.3125 0.0625]] 。。。。。。 Current round: 99 [[ 0.625 0.3125 0.0625]] Current round: 100 [[ 0.625 0.3125 0.0625]]
可以發現，從第60輪開始，我們的狀態概率分佈就不變了，一直保持在[0.625 0.3125 0.0625]，即62.5%的牛市，31.25%的熊市與6.25%的橫盤。那麼這個是巧合嗎？
我們現在換一個初始概率分佈試一試，現在我們用[0.7,0.1,0.2]作為初始概率分佈，然後這個狀態作為序列概率分佈的初始狀態$t_0$，將其帶入這個狀態轉移矩陣計算$t_1, t_2,t_3...$的狀態。程式碼如下：
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) vector1 = np.matrix([[0.7,0.1,0.2]], dtype=float) for i in range(100): vector1 = vector1*matrix print "Current round:" , i+1 print vector1
部分輸出結果如下：
Current round: 1 [[ 0.695 0.1825 0.1225]] Current round: 2 [[ 0.6835 0.22875 0.08775]] Current round: 3 [[ 0.6714 0.2562 0.0724]] Current round: 4 [[ 0.66079 0.273415 0.065795]] 。。。。。。。 Current round: 55 [[ 0.62500001 0.31249999 0.0625 ]] Current round: 56 [[ 0.62500001 0.31249999 0.0625 ]] Current round: 57 [[ 0.625 0.3125 0.0625]] 。。。。。。。 Current round: 99 [[ 0.625 0.3125 0.0625]] Current round: 100 [[ 0.625 0.3125 0.0625]]
可以看出，儘管這次我們採用了不同初始概率分佈，最終狀態的概率分佈趨於同一個穩定的概率分佈[0.625 0.3125 0.0625]， 也就是說我們的馬爾科夫鏈模型的狀態轉移矩陣收斂到的穩定概率分佈與我們的初始狀態概率分佈無關。這是一個非常好的性質，也就是說，如果我們得到了這個穩定概率分佈對應的馬爾科夫鏈模型的狀態轉移矩陣，則我們可以用任意的概率分佈樣本開始，帶入馬爾科夫鏈模型的狀態轉移矩陣，這樣經過一些序列的轉換，最終就可以得到符合對應穩定概率分佈的樣本。
這個性質不光對我們上面的狀態轉移矩陣有效，對於絕大多數的其他的馬爾科夫鏈模型的狀態轉移矩陣也有效。同時不光是離散狀態，連續狀態時也成立。
同時，對於一個確定的狀態轉移矩陣$P$，它的n次冪$P^n$在當n大於一定的值的時候也可以發現是確定的，我們還是以上面的例子為例，計算程式碼如下：
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) for i in range(10): matrix = matrix*matrix print "Current round:" , i+1 print matrix
輸出結果如下：
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]] 。。。。。。 Current round: 5 [[ 0.62502532 0.31247685 0.06249783] [ 0.6249537 0.31254233 0.06250397] [ 0.62497828 0.31251986 0.06250186]] Current round: 6 [[ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625]] Current round: 7 [[ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625]] 。。。。。。 Current round: 9 [[ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625]] Current round: 10 [[ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625] [ 0.625 0.3125 0.0625]]
我們可以發現，在$n \geq 6$以後，$P^n$的值穩定不再變化，而且每一行都為[0.625 0.3125 0.0625]，這和我們前面的穩定分佈是一致的。這個性質同樣不光是離散狀態，連續狀態時也成立。
好了，現在我們可以用數學語言總結下馬爾科夫鏈的收斂性質了：
如果一個非週期的馬爾科夫鏈有狀態轉移矩陣$P$, 並且它的任何兩個狀態是連通的，那麼$\lim_{n \to \infty}P_{ij}^n$與i無關，我們有：
1）$$\lim_{n \to \infty}P_{ij}^n = \pi(j)$$
2) $$\lim_{n \to \infty}P^n = \left( \begin{array}{ccc}
\pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \ldots&\ldots&\ldots&\ldots&\ldots \\ \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\
\ldots&\ldots&\ldots&\ldots&\ldots \end{array} \right)$$
3）$$\pi(j) = \sum\limits_{i=0}^{\infty}\pi(i)P_{ij}$$
4) $\pi$是方程$\pi P = \pi$的唯一非負解，其中：$$\pi = [\pi(1),\pi(2),...,\pi(j),...]\;\; \sum\limits_{i=0}^{\infty}\pi(i) = 1$$
上面的性質中需要解釋的有：
1）非週期的馬爾科夫鏈：這個主要是指馬爾科夫鏈的狀態轉化不是迴圈的，如果是迴圈的則永遠不會收斂。幸運的是我們遇到的馬爾科夫鏈一般都是非週期性的。用數學方式表述則是：對於任意某一狀態i，d為集合$\{n \mid n \geq 1,P_{ii}^n>0 \}$ 的最大公約數，如果 $d=1$ ，則該狀態為非週期的
2）任何兩個狀態是連通的：這個指的是從任意一個狀態可以通過有限步到達其他的任意一個狀態，不會出現條件概率一直為0導致不可達的情況。
3）馬爾科夫鏈的狀態數可以是有限的，也可以是無限的。因此可以用於連續概率分佈和離散概率分佈。、
4）$\pi$通常稱為馬爾科夫鏈的平穩分佈。
3. 基於馬爾科夫鏈取樣
如果我們得到了某個平穩分佈所對應的馬爾科夫鏈狀態轉移矩陣，我們就很容易採用出這個平穩分佈的樣本集。
假設我們任意初始的概率分佈是$\pi_0(x)$, 經過第一輪馬爾科夫鏈狀態轉移後的概率分佈是$\pi_1(x)$，。。。第i輪的概率分佈是$\pi_i(x)$。假設經過n輪後馬爾科夫鏈收斂到我們的平穩分佈$\pi(x)$，即：$$\pi_n(x) = \pi_{n+1}(x) = \pi_{n+2}(x) =... = \pi(x)$$
對於每個分佈$\pi_i(x)$，我們有：$$\pi_i(x) = \pi_{i-1}(x)P = \pi_{i-2}(x)P^2 = \pi_{0}(x)P^i$$
現在我們可以開始取樣了，首先，基於初始任意簡單概率分佈比如高斯分佈$\pi_0(x)$取樣得到狀態值$x_0$，基於條件概率分佈$P(x|x_0)$取樣狀態值$x_1$，一直進行下去，當狀態轉移進行到一定的次數時，比如到n次時，我們認為此時的取樣集$(x_n,x_{n+1},x_{n+2},...)$即是符合我們的平穩分佈的對應樣本集，可以用來做蒙特卡羅模擬求和了。
總結下基於馬爾科夫鏈的取樣過程：
1）輸入馬爾科夫鏈狀態轉移矩陣$P$，設定狀態轉移次數閾值$n_1$，需要的樣本個數$n_2$
2）從任意簡單概率分佈取樣得到初始狀態值$x_0$
3）for $t = 0$ to $n_1 + n_2 -1$: 從條件概率分佈$P(x|x_t)$中取樣得到樣本$x_{t+1}$
樣本集$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$即為我們需要的平穩分佈對應的樣本集。
4. 馬爾科夫鏈取樣小結
如果假定我們可以得到我們需要取樣樣本的平穩分佈所對應的馬爾科夫鏈狀態轉移矩陣，那麼我們就可以用馬爾科夫鏈取樣得到我們需要的樣本集，進而進行蒙特卡羅模擬。但是一個重要的問題是，隨意給定一個平穩分佈$\pi$,如何得到它所對應的馬爾科夫鏈狀態轉移矩陣$P$呢？這是個大問題。我們繞了一圈似乎還是沒有解決任意概率分佈取樣樣本集的問題。
幸運的是，MCMC取樣通過迂迴的方式解決了上面這個大問題，我們在下一篇來討論MCMC的取樣，以及它的使用改進版取樣: M-H取樣和Gibbs取樣.
（歡迎轉載，轉載請註明出處。歡迎溝通交流： liujianping-ok@163.com）