天天看點

python運作mcmc為何老出錯_Python可視化解析MCMC

馬爾可夫鍊可以定義為一個随機過程Y,其中t時刻各點的值隻取決于t-1時刻的值。這意味着随機過程在t時刻有狀态x的機率,給定它所有的過去狀态,等于在t時刻有狀态x的機率,給定它在t-1時刻的狀态。

python運作mcmc為何老出錯_Python可視化解析MCMC

如果可能的狀态集S是有限的,我們可以提供如下鍊的可視化表示:

python運作mcmc為何老出錯_Python可視化解析MCMC

每個圓圈表示一個狀态,在這種情況下,S={A, B, C},而箭頭表示我們的程序從一個狀态跳到另一個狀态的機率。我們可以把所有這些機率收集到一個矩陣P中,稱為過渡矩陣,如下:

python運作mcmc為何老出錯_Python可視化解析MCMC

在該情況下:

python運作mcmc為何老出錯_Python可視化解析MCMC

然後,在每個時間點上,我們可以描述(無條件)機率分布的過程, 它将是一個向量,其分量數等于S的維數。每個分量表示我們的過程取等于給定狀态的值的無條件機率。即:

python運作mcmc為何老出錯_Python可視化解析MCMC

關于μ的有趣性質是,它通過以下關系與過渡矩陣相連:

python運作mcmc為何老出錯_Python可視化解析MCMC

是以,一旦我們有一個已知的初始值的向量(這是有意義的,因為我們從一個可觀察到的狀态, 是以我們将有一個0的向量,在起始狀态的位置)可以計算在任何時候的分布過程。

此外,我們的向量有一個特定的值,這個等式成立:

python運作mcmc為何老出錯_Python可視化解析MCMC

如果存在這樣的值,我們調用相應的μ的不變分布的過程。

當我們讨論馬爾可夫鍊蒙特卡羅(MCMC)方法時,不變分布是一個關鍵的概念。後者包括一類從機率分布中抽樣的算法,它構造了一個以期望分布為不變分布的馬爾可夫鍊。

實際上,蒙特卡羅(MCMC)方法的目标是找到從不易取樣的分布中取樣的方法。為了繞過這個問題,有一些方法,比如拒絕抽樣和重要性抽樣,它們使用了一個更簡單的函數,叫做“proposal”。

讓我們模拟一個馬爾可夫鍊,考慮一個變量,其中今天的狀态可能隻取決于昨天的狀态。這個變量可能是天氣。讓我們考慮一下下面的鍊:

python運作mcmc為何老出錯_Python可視化解析MCMC

我們可以用和以前一樣的方法來解釋這個圖表。也就是說,如果今天是晴天,那麼明天還是晴天的機率是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運作mcmc為何老出錯_Python可視化解析MCMC

是以讓我們用Python實作它。為此,我設想檢查50天,然後我編碼Sunny = 1, Rainy = 2, Cloudy = 3。

m=np.zeros(150).reshape(50,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

Y[i] = 3

如果我繪制随機過程,我會得到這樣的結果:

python運作mcmc為何老出錯_Python可視化解析MCMC

有趣的是,如果我們計算機率分布清單的平均值(每個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])

我們從一個機率分布中建立了一個随機樣本它等于一個馬爾可夫鍊的不變分布。如果我們認為這個分布等于我們的目标分布(記住,很難從中取樣),我們就找到了一個繞過這個問題的方法。

原文連結:https://medium.com/analytics-vidhya/markov-chain-montecarlo-28dcde238e37