Python可视化解析MCMC

点击阅读原文，查看精彩日程！

作者 | Valentina Alto

译者 | 陆离

出品 | Python大本营（ID:pythonnews）

马尔可夫链（Markov Chain），又称为离散时间马尔可夫链，可以定义为一个随机过程Y，在某时间t上的任何一个点的值仅仅依赖于在时间t-1上的值。这就表示了我们的随机过程在时间t上具有状态x的概率，如果给出它之前所有的状态，那么就相当于在仅给出它在时间t-1的状态的时候，在时间t上具有状态x的概率。

如果可能的状态集S是有限的，那么，我们可以提供马尔可夫链的可视化表示结果，如下图所示：

上图中的每个圆圈都代表了一个状态，在这种情况下S={A, B, C}，而箭头则表示过程从一个状态跳到另一个状态的概率。我们可以在一个称为“转移矩阵”P中收集所有的这些概率数据，如下图所示：

那么，就有:

然后，在每个时间点上，我们可以描述过程的（无条件的）概率分布，这将是一个向量，其分量数等于S的维数。每个分量表示我们的过程取值等于给定状态的无条件概率。也就是:

关于上式中变量μ的比较有趣的性质是，它会通过以下等式的关系与转移矩阵相关联:

因此，一旦我们有了向量的已知初始值（这是可以理解的，因为我们是从一个可观察的状态开始的，因此将有一个包含多个0的向量，但在初始状态的位置上只有一个0），这样就可以计算过程在任何时间点上的分布了。

与此同时，我们的向量有一个特定的值，以使下面这个等式成立：

如果存在如上所述的一个值，我们将相应的变量μ称为过程的不变分布。

在讨论马尔可夫链蒙特卡罗（MCMC）方法的时候，不变分布是一个关键的概念。它包括一类从概率分布中抽样的算法，这个概率分布构造了一个马尔可夫链，而这个马尔可夫链则希望把这个分布作为它的不变分布。

事实上，蒙特卡罗方法的目标是要从不易抽样的分布中找到抽样的方法。要绕过这个问题，我们已有了一些方法，如拒绝抽样和重要性抽样等等，它们使用了一个更简单的函数，称为“proposal”（你可以点击这里找到相关的内容https://medium.com/datadriveninvestor/understanding-rejection-sampling-method-43b0006e0d0c）

让我们模拟一个马尔可夫链，现在，考虑一个变量，今天的状态可能只取决于昨天的状态，这个变量有可能指的是天气。所以让我们考虑下面的马尔可夫链：

我们可以用以前的方法来解释上图。也就是说，如果今天是晴天，则明天也是晴天的概率是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])

因此，我们从一个概率分布中创建了一个随机样本，而这个概率分布等于马尔可夫链的不变分布。如果我们认为这个分布等于目标分布（要记住，很难从中取样），那么就找到了绕过这个问题的办法。

原文地址

https://medium.com/analytics-vidhya/markov-chain-montecarlo-28dcde238e37

(本文由Python大本营翻译，转载请联系微信：1092722531)

◆

精彩推荐

◆

早鸟票限量5折！
2019 中国大数据技术大会（BDTC）再度来袭！豪华主席阵容及百位技术专家齐聚，15 场精选专题技术和行业论坛，超强干货+技术剖析+行业实践立体解读，深入解析热门技术在行业中的实践落地。

推荐阅读

你点的每个“在看”，我都认真当成了喜欢