通过给定的固定向量计算马尔可夫链
Compute Markov Chain by given stationary vector
我正在做 Monte Carlo 课程的作业,我被要求找到一个具有 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!
算法方法
计算马尔可夫链平稳分布的标准方法是求解线性方程,例如如此处所述 https://stephens999.github.io/fiveMinuteStats/stationary_distribution.html。
您的问题的解决方案是相反的 - 解决相同的方程式,除了在您的情况下,您具有平稳分布但没有过渡 probabilities/rates。
但是,这种方法的问题是您可能会构造一个线性方程组,其中的变量比方程多得多。这严重减少了您对构建的马尔可夫链的拓扑结构的选择。
幸运的是,您似乎对构造的马尔可夫链的拓扑结构没有任何限制,因此您可以做出一些妥协。你可以做的是禁用大多数转换,即给它们零 probability/rate,并且每个状态只启用一个转换。这可能会产生某种环形拓扑,但它应该确保您的线性方程组有解。
原始示例
考虑平稳分布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 然后回到 x 概率为 1。
请注意,拟合解的数量可能是无限的(即使对于示例中所示的简化线性方程组也是如此),例如在这种情况下,probabilities/rates 可以是任何正值,只要它们都相等即可。
所以,一步步解决
- 按照描述构建线性方程组。
- 求解构造的线性方程组
- 您构造的线性方程组的解描述了您正在寻找的马尔可夫链。如果需要,可以简单地重建整个转换矩阵。
我正在做 Monte Carlo 课程的作业,我被要求找到一个具有 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!
算法方法
计算马尔可夫链平稳分布的标准方法是求解线性方程,例如如此处所述 https://stephens999.github.io/fiveMinuteStats/stationary_distribution.html。 您的问题的解决方案是相反的 - 解决相同的方程式,除了在您的情况下,您具有平稳分布但没有过渡 probabilities/rates。
但是,这种方法的问题是您可能会构造一个线性方程组,其中的变量比方程多得多。这严重减少了您对构建的马尔可夫链的拓扑结构的选择。 幸运的是,您似乎对构造的马尔可夫链的拓扑结构没有任何限制,因此您可以做出一些妥协。你可以做的是禁用大多数转换,即给它们零 probability/rate,并且每个状态只启用一个转换。这可能会产生某种环形拓扑,但它应该确保您的线性方程组有解。
原始示例
考虑平稳分布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 然后回到 x 概率为 1。 请注意,拟合解的数量可能是无限的(即使对于示例中所示的简化线性方程组也是如此),例如在这种情况下,probabilities/rates 可以是任何正值,只要它们都相等即可。
所以,一步步解决
- 按照描述构建线性方程组。
- 求解构造的线性方程组
- 您构造的线性方程组的解描述了您正在寻找的马尔可夫链。如果需要,可以简单地重建整个转换矩阵。