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

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

如果可能的狀態(tài)集S是有限的，那么，我們可以提供馬爾可夫鏈的可視化表示結(jié)果，如下圖所示：

上圖中的每個圓圈都代表了一個狀態(tài)，在這種情況下S={A, B, C}，而箭頭則表示過程從一個狀態(tài)跳到另一個狀態(tài)的概率。我們可以在一個稱為“轉(zhuǎn)移矩陣”P中收集所有的這些概率數(shù)據(jù)，如下圖所示：

那么，就有:

然后，在每個時間點上，我們可以描述過程的（無條件的）概率分布，這將是一個向量，其分量數(shù)等于S的維數(shù)。每個分量表示我們的過程取值等于給定狀態(tài)的無條件概率。也就是:

關(guān)于上式中變量μ的比較有趣的性質(zhì)是，它會通過以下等式的關(guān)系與轉(zhuǎn)移矩陣相關(guān)聯(lián):

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

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

如果存在如上所述的一個值，我們將相應(yīng)的變量μ稱為過程的不變分布。

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

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

讓我們模擬一個馬爾可夫鏈，現(xiàn)在，考慮一個變量，今天的狀態(tài)可能只取決于昨天的狀態(tài)，這個變量有可能指的是天氣。所以讓我們考慮下面的馬爾可夫鏈：

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

我們可以在以下的轉(zhuǎn)移矩陣中收集表示上圖中箭頭的數(shù)組:

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 ]])

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

由于我們有一個初始的變量μ和一個轉(zhuǎn)移矩陣，因此就可以在任意時間點t上計算μ的值。因此，有了這些之后，我想根據(jù)每個t值的概率分布來創(chuàng)建一個隨機(jī)過程（具有馬爾可夫鏈的屬性，因此可以只依賴于前一個時間段）。

這意味著我得到的隨機(jī)變量Y將會有一些等于瞬間數(shù)量的分量，而每個分量都是根據(jù)瞬間的概率分布來實現(xiàn)的過程。為此，我們希望從均勻分布中生成一個隨機(jī)數(shù)，并設(shè)置如下規(guī)則：

讓我們用Python語言來實現(xiàn)程序代碼。為此，我假設(shè)了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

如果我用圖表來繪制隨機(jī)過程，將會得到以下類似的結(jié)果：

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

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

這近似于不變分布，它可以進(jìn)行如下的計算:

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])

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