python运行mcmc为何老出错_Python可视化解析MCMC

马尔可夫链可以定义为一个随机过程Y，其中t时刻各点的值只取决于t-1时刻的值。这意味着随机过程在t时刻有状态x的概率，给定它所有的过去状态，等于在t时刻有状态x的概率，给定它在t-1时刻的状态。

如果可能的状态集S是有限的，我们可以提供如下链的可视化表示:

每个圆圈表示一个状态，在这种情况下，S={A, B, C}，而箭头表示我们的进程从一个状态跳到另一个状态的概率。我们可以把所有这些概率收集到一个矩阵P中，称为过渡矩阵，如下:

在该情况下：

然后,在每个时间点上,我们可以描述(无条件)概率分布的过程, 它将是一个向量，其分量数等于S的维数。每个分量表示我们的过程取等于给定状态的值的无条件概率。即：

关于μ的有趣性质是，它通过以下关系与过渡矩阵相连：

所以,一旦我们有一个已知的初始值的向量(这是有意义的,因为我们从一个可观察到的状态, 因此我们将有一个0的向量，在起始状态的位置)可以计算在任何时候的分布过程。

此外，我们的向量有一个特定的值，这个等式成立:

如果存在这样的值,我们调用相应的μ的不变分布的过程。

当我们讨论马尔可夫链蒙特卡罗(MCMC)方法时，不变分布是一个关键的概念。后者包括一类从概率分布中抽样的算法，它构造了一个以期望分布为不变分布的马尔可夫链。

实际上，蒙特卡罗(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=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

如果我绘制随机过程，我会得到这样的结果：

有趣的是，如果我们计算概率分布列表的平均值(每个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