通过给定的平稳向量计算马尔可夫链

Question

我正在蒙特卡洛课程做作业，我被要求找到一个具有 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 中！

Answer 1

算法方法

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

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

原始示例

考虑平稳分布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。 请注意，拟合解的数量可能是无限的（即使对于示例中所示的简化线性方程组也是如此），例如在这种情况下，概率/比率可以是任何正值，只要它们都相等。

所以，一步一步解决

1. 按照描述构建线性方程组。
2. 求解构造的线性方程组
3. 您构建的线性方程组的解描述了您正在寻找的马尔可夫链。如果您愿意，可以简单地重建整个转换矩阵。