问题描述
我正在尝试计算一个计算马尔可夫链均衡概率的函数。 对于这个问题，已经得到了我的转移矩阵。
现在我正在尝试定义一个名为 Jacobi 的函数，但对最有效的方法感到困惑。有关如何执行此操作的任何建议？
到目前为止，我已尝试将其设置为方程组并求解 x=a^(-1)*b，但由于转换矩阵是奇异的，因此无法正确实现。
我知道我需要将转换矩阵乘以一个变量矩阵才能得到 7 个单独的方程。然后我需要添加方程 x0 + x1 + x2 + x3 + x4 + x5 + x6 = 1。在我得到所有 8 个方程后，我可以求解 x0 到 x6 以获得我的均衡概率。你知道我如何在 python 代码中实现这个过程吗？
解决方法
我不确定雅可比方法是否适用于马尔可夫转换矩阵。不过，还有其他更简单的技术可以找到平稳分布。首先像你描述的那样求解方程组：
import numpy as np
>>> M = np.array([[0.084,0.1008,0.168,0.224,0.1494,0.0498],[0.084,[0.5768,0.0498,0.0,0.0],[0.3528,[0.1848,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.02876204])