查看原文
其他

Python可视化解析MCMC

Valentina Alto Python大本营 2019-10-30

点击阅读原文,查看精彩日程!


作者 | 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.337777780.293333330.36888889])


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


原文地址

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


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


精彩推荐


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


推荐阅读


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

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存