专栏首页谓之小一MCMC之马尔可夫链

MCMC之马尔可夫链

MCMC之蒙特卡罗方法之

1.马尔可夫链简介

因为某一时刻状态转移只依赖于它的前一个状态,那么我们只要能求出系统中任意两个状态之间的转移概率,进而得到状态转移概率矩阵,那么马尔科夫链的模型便定了。以下图股市模型为例,共有三个状态,分别为牛市(Bull market)、熊市(Bear market)、横盘(Stagnant market)。每一个状态都能够以一定概率转移到下一状态,比如牛市以0.075的概率转移到横盘的概率,这些状态转移概率图可以转换为矩阵的形式进行表示。

如果我们定义矩阵$P$某一位置P(i,j)的值为P(j|i),即从状态i转移到状态j的概率,并定义牛市的状态为0、熊市状态为1、横盘状态为2,这样便得到马尔可夫链模型的状态转移矩阵。

那么马尔科夫链模型的状态转移矩阵和蒙特卡罗方法所需要的概率分布样本集有什么关系呢?

2.马尔可夫链状态转移矩阵性质

得到马尔可夫链状态转移矩阵,我们看看马尔可夫链模型状态转移矩阵的性质。仍以上面的状态转移矩阵为例,假设当前股市的概率分布为[0.3, 0.4, 0.3],即30%概率的牛市、40%概率的熊市、30%概率的横盘。如果以这个状态作为序列概率分布的初始状态t0,与状态转移矩阵计算得到t1,t2,t3,…时刻的概率,相应代码和结果如下。

import numpy as np


def markov_chain():
    matrix = np.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05],
                        [0.25, 0.25, 0.5]], dtype=float)
    current = np.matrix([[0.3, 0.4, 0.3]], dtype=float)
    for i in range(100):
        current = current * matrix
        print "Current round:", i + 1
        print current


if __name__ == '__main__':
    markov_chain()
Current round: 1
[[0.405  0.4175 0.1775]]
Current round: 2
[[0.4715  0.40875 0.11975]]
Current round: 3
[[0.5156 0.3923 0.0921]]
Current round: 4
[[0.54591  0.375535 0.078555]]
...
...
Current round: 58
[[0.62499999 0.31250001 0.0625    ]]
Current round: 59
[[0.62499999 0.3125     0.0625    ]]
Current round: 60
[[0.625  0.3125 0.0625]]
Current round: 61
[[0.625  0.3125 0.0625]]
Current round: 62
[[0.625  0.3125 0.0625]]
...
...
Current round: 98
[[0.625  0.3125 0.0625]]
Current round: 99
[[0.625  0.3125 0.0625]]
Current round: 100
[[0.625  0.3125 0.0625]]

可以发现从第60轮开始,状态转移矩阵的概率分布就不变了,一直保持在[0.625, 0.3125, 0.0625],即62.5%的概率的牛市、31.25%概率的熊市、6.25%概率的横盘,那么这个是巧合吗?

答案并不是巧合,如果我们换一个初始概率分布t0,比如以[0.7, 0.1, 0.2]作为初始概率分布,然后带入状态转移矩阵得到t1,t2,t3,…时刻的概率,会发现得到同样的结果。即最终的状态概率分布会趋于同一个稳定概率分布[0.625, 0.3125, 0.0625],也就是说,马尔可夫链的状态转移矩阵收敛到稳定概率分布与初始状态概率分布无关。

上述结果是一个非常好的形式,比如我们得到了稳定概率分布所对应的马尔可夫链模型的状态转移矩阵,那么可以用任意的概率分布样本开始,带入马尔可夫链状态转移矩阵,然后就可以得到符合对应稳定概率分布的样本。这个性质不光对于上面的状态转移矩阵有效,对于绝大多数的其他马尔可夫链模型的状态转移矩阵也有效。同时不光是离散状态,连续状态情况下也成立。

import numpy as np


def markov_chain():
    matrix = np.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]], dtype=float)
    for i in range(10):
        matrix = matrix * matrix
        print "Current round:", i + 1
        print matrix


if __name__ == '__main__':
    markov_chain()
Current round: 1
[[0.8275  0.13375 0.03875]
 [0.2675  0.66375 0.06875]
 [0.3875  0.34375 0.26875]]
Current round: 2
[[0.73555  0.212775 0.051675]
 [0.42555  0.499975 0.074475]
 [0.51675  0.372375 0.110875]]
...
Current round: 5
[[0.62502532 0.31247685 0.06249783]
 [0.6249537  0.31254233 0.06250397]
 [0.62497828 0.31251986 0.06250186]]
Current round: 6
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]
...
Current round: 10
[[0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]
 [0.625  0.3125 0.0625]]

上面性质中,有几点需要特意说明

3.基于马尔可夫链采样

4.马尔可夫链总结

如果假定我们可以得到所需要采样样本的平稳分布所对应的马尔可夫链状态转移矩阵,那么我们就可以用马尔可夫链采样得到我们需要的样本集,进而进行蒙特卡罗模拟。但是现在还有个很重要的问题,随意给定一个平稳分布π ,如何得到它所对应的马尔可夫链状态转移矩阵P呢?下篇文章,我们将重点介绍MCMC采用通过与会的方式解决上述问题,以及改进版的M-H采样和Gibbs采样。

你看到的这篇文章来自于公众号「谓之小一」,欢迎关注我阅读更多文章。

本文分享自微信公众号 - 谓之小一(weizhixiaoyi),作者:谓之小一

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-12-16

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 电影知识图谱问答(四)| 问句理解及答案推理

    上篇文章《电影知识图谱问答(三)|Apache Jena知识存储及SPARQL知识检索》中讲到如何将处理后的RDF数据存储至Apache Jena数据库之中、如...

    小一
  • 机器学习之随机森林

    随机森林(Random Forest)是一个非常灵活的机器学习方法,从市场营销到医疗保险有着众多的应用。例如用于市场营销对客户获取和存留建模或预测病人的疾病风险...

    小一
  • 机器学习之分类与回归树(CART)

    分类与回归树的英文是Classfication And Regression Tree,缩写为CART。CART算法采用二分递归分割的技术将当前样本集分为两个子...

    小一
  • AI未来发展的三种模式 | CCF-GAIR人工智能前沿专场随笔

    2017年7月7日-2017年7月9日,由中国计算机协会CCF主办、雷锋网和香港中文大学(深圳)联合承办的CCF-GAIR全球人工智能与机器人峰会,在深圳大中华...

    数据派THU
  • 告诉你做数据分析必须学R的4个理由

    论坛君:你很可能已经听说过 R,或许你知道 R 是一种编程语言,而且知道它与统计学有关,但它是否适合您呢?本文作者将试图向大家讲解他对R的看法,分享他认为试用开...

    小莹莹
  • 如何批量下载电视剧

    生活中经常会有下载点东西的需要,有些网页虽然有批量下载的功能,但很多时候都不好用。并且我觉得网站可能就想让你多点几下,在网站上多停留一会。于其在网站上点来点去,...

    数据处理与分析
  • 界面心理学是怎么一回事?

    汽车仪表盘设计的演变 ? Dashboard 起源于汽车的仪表显示盘,在国外通常被称为仪表盘,国内从一开始的“后台”到现在称之为“信息面板”,得之于web后台到...

    BestSDK
  • 玩转JS的类型转换黑科技0.前言1.奇葩例子2.从[]==![]开始3.从已有的得到想不到的4.关于(a==1 && a==2 && a==3)4.2 ===

    js身为一种弱类型的语言,不用像c语言那样要定义int、float、double、string等等数据类型,因为允许变量类型的隐式转换和允许强制类型转换。我们在...

    lhyt
  • Vue-cli 目录结构

    hedeqiang
  • Redis(8)——发布/订阅与Stream

    发布/ 订阅系统 是 Web 系统中比较常用的一个功能。简单点说就是 发布者发布消息,订阅者接受消息,这有点类似于我们的报纸/ 杂志社之类的: (借用前边的一张...

    我没有三颗心脏

扫码关注云+社区

领取腾讯云代金券