龙空技术网

马尔可夫链你知道多少?Python可视化解析MCMC

千锋IT教育 1229

前言:

如今姐妹们对“mc算法 蒙特卡罗算法”大概比较珍视,朋友们都需要知道一些“mc算法 蒙特卡罗算法”的相关文章。那么小编也在网上收集了一些关于“mc算法 蒙特卡罗算法””的相关内容,希望咱们能喜欢,同学们一起来学习一下吧!

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

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

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

那么,就有:

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

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

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

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

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

在讨论马尔可夫链蒙特卡罗(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[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])

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

标签: #mc算法 蒙特卡罗算法 #python mcmc