馬爾可夫鏈你知道多少？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])

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