馬爾可夫鏈（Markov Chain），又稱為離散時(shí)間馬爾可夫鏈，可以定義為一個(gè)隨機(jī)過(guò)程Y，在某時(shí)間t上的任何一個(gè)點(diǎn)的值僅僅依賴于在時(shí)間t-1上的值。這就表示了我們的隨機(jī)過(guò)程在時(shí)間t上具有狀態(tài)x的概率，如果給出它之前所有的狀態(tài)，那么就相當(dāng)于在僅給出它在時(shí)間t-1的狀態(tài)的時(shí)候，在時(shí)間t上具有狀態(tài)x的概率。
如果可能的狀態(tài)集S是有限的，那么，我們可以提供馬爾可夫鏈的可視化表示結(jié)果，如下圖所示：
上圖中的每個(gè)圓圈都代表了一個(gè)狀態(tài)，在這種情況下S={A, B, C}，而箭頭則表示過(guò)程從一個(gè)狀態(tài)跳到另一個(gè)狀態(tài)的概率。我們可以在一個(gè)稱為“轉(zhuǎn)移矩陣”P中收集所有的這些概率數(shù)據(jù)，如下圖所示：
那么，就有:
然后，在每個(gè)時(shí)間點(diǎn)上，我們可以描述過(guò)程的（無(wú)條件的）概率分布，這將是一個(gè)向量，其分量數(shù)等于S的維數(shù)。每個(gè)分量表示我們的過(guò)程取值等于給定狀態(tài)的無(wú)條件概率。也就是:
關(guān)于上式中變量μ的比較有趣的性質(zhì)是，它會(huì)通過(guò)以下等式的關(guān)系與轉(zhuǎn)移矩陣相關(guān)聯(lián):
因此，一旦我們有了向量的已知初始值（這是可以理解的，因?yàn)槲覀兪菑囊粋€(gè)可觀察的狀態(tài)開始的，因此將有一個(gè)包含多個(gè)0的向量，但在初始狀態(tài)的位置上只有一個(gè)0），這樣就可以計(jì)算過(guò)程在任何時(shí)間點(diǎn)上的分布了。
與此同時(shí)，我們的向量有一個(gè)特定的值，以使下面這個(gè)等式成立：
如果存在如上所述的一個(gè)值，我們將相應(yīng)的變量μ稱為過(guò)程的不變分布。
在討論馬爾可夫鏈蒙特卡羅（MCMC）方法的時(shí)候，不變分布是一個(gè)關(guān)鍵的概念。它包括一類從概率分布中抽樣的算法，這個(gè)概率分布構(gòu)造了一個(gè)馬爾可夫鏈，而這個(gè)馬爾可夫鏈則希望把這個(gè)分布作為它的不變分布。
事實(shí)上，蒙特卡羅方法的目標(biāo)是要從不易抽樣的分布中找到抽樣的方法。要繞過(guò)這個(gè)問(wèn)題，我們已有了一些方法，如拒絕抽樣和重要性抽樣等等，它們使用了一個(gè)更簡(jiǎn)單的函數(shù)，稱為“proposal”
讓我們模擬一個(gè)馬爾可夫鏈，現(xiàn)在，考慮一個(gè)變量，今天的狀態(tài)可能只取決于昨天的狀態(tài)，這個(gè)變量有可能指的是天氣。所以讓我們考慮下面的馬爾可夫鏈：
我們可以用以前的方法來(lái)解釋上圖。也就是說(shuō)，如果今天是晴天，則明天也是晴天的概率是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 ]])
另外，也有一個(gè)初始值，比如說(shuō)“多云”，因此我們已經(jīng)有了y的初始分布，即μ _0=[0,0,1]。
由于我們有一個(gè)初始的變量μ和一個(gè)轉(zhuǎn)移矩陣，因此就可以在任意時(shí)間點(diǎn)t上計(jì)算μ的值。因此，有了這些之后，我想根據(jù)每個(gè)t值的概率分布來(lái)創(chuàng)建一個(gè)隨機(jī)過(guò)程（具有馬爾可夫鏈的屬性，因此可以只依賴于前一個(gè)時(shí)間段）。
這意味著我得到的隨機(jī)變量Y將會(huì)有一些等于瞬間數(shù)量的分量，而每個(gè)分量都是根據(jù)瞬間的概率分布來(lái)實(shí)現(xiàn)的過(guò)程。為此，我們希望從均勻分布中生成一個(gè)隨機(jī)數(shù)，并設(shè)置如下規(guī)則：
讓我們用Python語(yǔ)言來(lái)實(shí)現(xiàn)程序代碼。為此，我假設(shè)了50天的測(cè)試，然后我輸入：
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
如果我用圖表來(lái)繪制隨機(jī)過(guò)程，將會(huì)得到以下類似的結(jié)果：
在這個(gè)過(guò)程中比較有趣的是，如果計(jì)算這些概率分布中列表的平均值（每個(gè)t值對(duì)應(yīng)一個(gè)），我們將會(huì)得到：
[np.mean(m[:,0]), np.mean(m[:,1]), np.mean(m[:,2])] Output: [0.3239190123456788, 0.2888770370370369, 0.3872039506172838]
這近似于不變分布，它可以進(jìn)行如下的計(jì)算:
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])
因此，我們從一個(gè)概率分布中創(chuàng)建了一個(gè)隨機(jī)樣本，而這個(gè)概率分布等于馬爾可夫鏈的不變分布。如果我們認(rèn)為這個(gè)分布等于目標(biāo)分布（要記住，很難從中取樣），那么就找到了繞過(guò)這個(gè)問(wèn)題的辦法。
聯(lián)系客服
