馬爾可夫鏈你知道多少？Python視覺化解析MCMC

馬爾可夫鏈（Markov Chain），又稱為離散時間馬爾可夫鏈，可以定義為一個隨機過程Y，在某時間t上的任何一個點的值僅僅依賴於在時間t-1上的值。這就表示了我們的隨機過程在時間t上具有狀態x的機率，如果給出它之前所有的狀態，那麼就相當於在僅給出它在時間t-1的狀態的時候，在時間t上具有狀態x的機率。

如果可能的狀態集S是有限的，那麼，我們可以提供馬爾可夫鏈的視覺化表示結果，如下圖所示：

上圖中的每個圓圈都代表了一個狀態，在這種情況下S={A, B, C}，而箭頭則表示過程從一個狀態跳到另一個狀態的機率。我們可以在一個稱為“轉移矩陣”P中收集所有的這些機率資料，如下圖所示：

那麼，就有:

然後，在每個時間點上，我們可以描述過程的（無條件的）機率分佈，這將是一個向量，其分量數等於S的維數。每個分量表示我們的過程取值等於給定狀態的無條件機率。也就是:

關於上式中變數μ的比較有趣的性質是，它會透過以下等式的關係與轉移矩陣相關聯:

因此，一旦我們有了向量的已知初始值（這是可以理解的，因為我們是從一個可觀察的狀態開始的，因此將有一個包含多個0的向量，但在初始狀態的位置上只有一個0），這樣就可以計算過程在任何時間點上的分佈了。

與此同時，我們的向量有一個特定的值，以使下面這個等式成立：

如果存在如上所述的一個值，我們將相應的變數μ稱為過程的不變分佈。

在討論馬爾可夫鏈蒙特卡羅（MCMC）方法的時候，不變分佈是一個關鍵的概念。它包括一類從機率分佈中抽樣的演算法，這個機率分佈構造了一個馬爾可夫鏈，而這個馬爾可夫鏈則希望把這個分佈作為它的不變分佈。

事實上，蒙特卡羅方法的目標是要從不易抽樣的分佈中找到抽樣的方法。要繞過這個問題，我們已有了一些方法，如拒絕抽樣和重要性抽樣等等，它們使用了一個更簡單的函式，稱為“proposal”

讓我們模擬一個馬爾可夫鏈，現在，考慮一個變數，今天的狀態可能只取決於昨天的狀態，這個變數有可能指的是天氣。所以讓我們考慮下面的馬爾可夫鏈：

我們可以用以前的方法來解釋上圖。也就是說，如果今天是晴天，則明天也是晴天的機率是50%，而下雨的機率是15%，是多雲天氣的機率是35%。

我們可以在以下的轉移矩陣中收集表示上圖中箭頭的陣列:

import numpy as np P = np.array([[0.5, 0.15, 0.35], [0.45, 0.45, 0.1], [0.1, 0.3, 0.6]]) P Output: array([[0.5 , 0.15, 0.35], [0.45, 0.45, 0.1 ], [0.1 , 0.3 , 0.6 ]])

另外，也有一個初始值，比如說“多雲”，因此我們已經有了y的初始分佈，即μ _0=[0,0,1]。

由於我們有一個初始的變數μ和一個轉移矩陣，因此就可以在任意時間點t上計算μ的值。因此，有了這些之後，我想根據每個t值的機率分佈來建立一個隨機過程（具有馬爾可夫鏈的屬性，因此可以只依賴於前一個時間段）。

這意味著我得到的隨機變數Y將會有一些等於瞬間數量的分量，而每個分量都是根據瞬間的機率分佈來實現的過程。為此，我們希望從均勻分佈中生成一個隨機數，並設定如下規則：

讓我們用Python語言來實現程式程式碼。為此，我假設了50天的測試，然後我輸入：

Sunny = 1, Rainy = 2, Cloudy = 3.
m=np.zeros(150).reshape(50,3) m[0]=[0,0,1] ndays = 50 Y=[0]*ndays u = np.random.uniform(0,1,50) for i in range(1, ndays): tmp=[] m[i] = m[i-1].dot(P) if u[i] < m[i][0]: Y[i]=1 elif u[i] < m[i][0] + m[i][1]: Y[i] = 2 else: Y[i] = 3

如果我用圖表來繪製隨機過程，將會得到以下類似的結果：

在這個過程中比較有趣的是，如果計算這些機率分佈中列表的平均值（每個t值對應一個），我們將會得到：

[np.mean(m[:,0]), np.mean(m[:,1]), np.mean(m[:,2])] Output: [0.3239190123456788, 0.2888770370370369, 0.3872039506172838]

這近似於不變分佈，它可以進行如下的計算:

a=np.array([[-0.5, 0.45, 0.1], [0.15, -0.55, 0.3], [1,1,1]]) b=np.array([0,0,1]) mu = np.linalg.solve(a, b) mu Output: array([0.33777778, 0.29333333, 0.36888889])

因此，我們從一個機率分佈中建立了一個隨機樣本，而這個機率分佈等於馬爾可夫鏈的不變分佈。如果我們認為這個分佈等於目標分佈（要記住，很難從中取樣），那麼就找到了繞過這個問題的辦法。