python 可视化_Python可视化解析MCMC

作者 | Valentina Alto 译者 | 陆离 出品 | Python大本营(ID:pythonnews)

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

24bb912a08c2634a157891915611f63c.png

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

21a797f6937ace50fe2ee8b4f319190c.png

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

40123ffefabedaeac256de70969efbff.png

那么,就有:

d7c2dd03172d31f4889c41965957175b.png

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

9d051a81766a9f02d42d2453e4dd4744.png

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

e8eb6ef029c1b6f527bd57c9f3238ef0.png

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

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

8fd9ccd5ccc26c788579e723502746c1.png

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

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

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

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

83053ac6f3b828c6a187c30f046373d8.png

我们可以用以前的方法来解释上图。也就是说,如果今天是晴天,则明天也是晴天的概率是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将会有一些等于瞬间数量的分量,而每个分量都是根据瞬间的概率分布来实现的过程。为此,我们希望从均匀分布中生成一个随机数,并设置如下规则:

02166249683aca5d3e94fbcf6d794692.png

让我们用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

如果我用图表来绘制随机过程,将会得到以下类似的结果:

5b572a40fcd6e66acda8aff39d8e47dc.png

在这个过程中比较有趣的是,如果计算这些概率分布中列表的平均值(每个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 场 精选 专题技术和行业论坛, 超强干货+技术剖析+行业实践立体解读 ,深入解析热门技术在行业中的实践落地。

72b821235a7a60a515a4ce66ed3a83c5.png

推荐阅读

  • 每个新手程序员都应该知道的Python开发技巧

  • Python带你分析,英超是否已经大结局
  • Python GUI开发,效率提升10倍的方法!
  • Python老司机给上路新手的3点忠告

  • Python新工具:用三行代码提取PDF表格数据
  • 微软也爱Python!VS Code Python全新发布!Jupyter Notebook原生支持终于来了!
  • 5大必知的图算法,附Python代码实现

  • 吐血整理!140种Python标准库、第三方库和外部工具都有了

  • 如何用爬虫技术帮助孩子秒到心仪的幼儿园(基础篇)

  • Python传奇:30年崛起之路

6cabe13d40e050eb0b90120ffecc47c2.png 你点的每个“在看”,我都认真当成了喜欢
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值