-1

我正在尝试计算一个计算马尔可夫链平衡概率的函数。对于这个问题，我已经得到了我的转换矩阵。

现在我正在尝试定义一个名为 Jacobi 的函数，但对最有效的方法感到困惑。关于如何做到这一点的任何建议？

到目前为止，我已经尝试将它设置为方程组并求解 x=a^(-1)*b 但由于转换矩阵是奇异的而无法正确实现它。

我知道我需要将转换矩阵乘以一个变量矩阵来获得 7 个单独的方程。然后我需要添加方程 x0 + x1 + x2 + x3 + x4 + x5 + x6 = 1。在我拥有所有 8 个方程后，我可以求解 x0 到 x6 以获得我的平衡概率。你知道我如何在 python 代码中实现这个过程吗？

4

1 回答 1

1

我不确定 Jacobi 方法是否适用于马尔可夫转移矩阵。不过，还有其他更简单的技术可以找到平稳分布。首先通过求解您所描述的方程组：

import numpy as np


>>> M = np.array([[0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.5768, 0.224, 0.1494, 0.0498, 0.0, 0.0, 0.0], 
                  [0.3528, 0.224, 0.224, 0.1494, 0.0498, 0.0, 0.0], 
                  [0.1848, 0.168, 0.224, 0.224, 0.1494, 0.0498, 0.0], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498]])
>>>
>>> I = np.eye(len(M)) # identity matrix
>>>
>>> A = M.T - I  # rearrange each "equation" in the system
>>>
>>> A = np.vstack([A, np.ones((1, len(A)))]) # add row of ones
>>>
>>> b = np.zeros(len(A)) # prepare solution vector
>>> b[-1] = 1.0 # add a one to the last entry
>>>
>>> np.linalg.solve(A.T @ A, A.T @ b) # solve the normal equation
array([0.22288974, 0.14776044, 0.1781388 , 0.181211  , 0.1504296 ,
       0.09080838, 0.02876204])

您也可以反复乘以M自身，直到它收敛：

>>> np.linalg.matrix_power(M, 25)[0]
array([0.22288974, 0.14776044, 0.1781388 , 0.181211  , 0.1504296 ,
       0.09080838, 0.02876204])
于 2021-04-20T05:31:55.627 回答