用给定的静止向量计算马尔可夫链

问题描述 投票:0回答:1

我正在蒙特卡洛课程中做作业,我被要求找到一个马尔可夫链的矩阵,有6个状态,分别为0,1,2,3,4,5,这样经过足够长的一段时间后,我们花费的时间与在每个州中的数字5,10,5,10,25,60。

我知道如果我们有转换矩阵,这就是我们得到的静止向量。我必须使用Metropolis算法,但我发现的所有解释和示例都基于Metropolis-Hasting算法。

我有的算法的伪代码:

select x
 Loop over repetitions t=1,2...
 select y from Nx using density Gx
 put h=min(1, f(y)/f(x))
 if U ~ U(0, 1) < h then x <- y
end loop

我正在寻找一步一步解释如何为给定问题实现算法,最好是在python中!

python markov-chains mcmc
1个回答
1
投票

算法方法

计算马尔可夫链的平稳分布的标准方法是线性方程的解,例如,如此处所述https://stephens999.github.io/fiveMinuteStats/stationary_distribution.html。你的问题的解决方案是相反的 - 解决相同的方程,除了在你的情况下,你有固定的分布,但你没有转换概率/率。

然而,这种方法的问题在于,您可以构建一个线性方程组,其变量比方程式多得多。这严重降低了您对构建马尔可夫链拓扑的选择。幸运的是,您似乎对构造的马尔可夫链的拓扑结构没有任何限制,因此您可以做出一些妥协。你可以做的是禁用大多数转换,即给它们零概率/速率,并且每个状态只启用一个转换。这可能会产生某种环形拓扑,但它应该确保您的线性方程组有一个解决方案。

Primitive example

考虑平稳分布Pi =(x = 1/3,y = 1/3,z = 1/3)

构造你的线性方程组

Pi(x) = 1/3 =  Pr(y,x) * Pi(y)
Pi(y) = 1/3 =  Pr(z,y) * Pi(z)
Pi(z) = 1/3 =  Pr(x,z) * Pi(x)

在这种情况下,解是Pr(y,x)= Pr(z,y)= Pr(x,z)= 1并且所获得的马尔可夫链刚刚从x到z循环到y并且以概率1回到x。注意,拟合解的数量可以是无限的(即使对于如示例中所示的简化的线性方程组),例如在这种情况下,概率/速率可以是任何正值,只要它们都相等即可。

So, step by step solution

  1. 构建如所描述的线性方程组。
  2. 求解构造的线性方程组
  3. 您构建的线性方程组的解决方案描述了您正在寻找的马尔可夫链。如果您愿意,可以轻松地重建整个转换矩阵。
© www.soinside.com 2019 - 2024. All rights reserved.