概率矩阵 Numpy 的终端概率

Question

我有一个矩阵 m，它表示从状态转换到状态的概率。

例如对于下面的示例，我将始终卡在状态 1、3、4，而状态 2 我将随机转换到 4 个状态之一。

import numpy as np
m = np.eye(4)
m[1] = 0.25
print(m)

[[1.   0.   0.   0.  ]

 [0.25 0.25 0.25 0.25]

 [0.   0.   1.   0.  ]

 [0.   0.   0.   1.  ]]

如何找到表示无限转换后最终状态的矩阵？

例如如果我这样做，我会得到状态 1,3,4 --> 100% 坚持 1,3,4 但状态 2 --> 1/3 机会最终出现在所有其他状态的直观结果。由于来自状态 2 的所有案例最终通过多次转换在 1、3、4 之间平均分配。

t = m
for _ in range(100_000):
    t = t @ t
print(t)

[[1.         0.         0.         0.        ]

 [0.33333333 0.         0.33333333 0.33333333]

 [0.         0.         1.         0.        ]

 [0.         0.         0.         1.        ]]

如何在不使用重复乘法的情况下进行计算？我认为它对应于矩阵的特征向量/特征值，但是当我计算它时，我得到了一些非常不同的东西。

np.linalg.eig(m)

[[0.        , 0.9486833 , 0.        , 0.        ],

[1.        , 0.31622777, 0.31622777, 0.31622777],

[0.        , 0.        , 0.9486833 , 0.        ],

[0.        , 0.        , 0.        , 0.9486833 ]]

有没有一种方法可以使用 numpy 进行计算？我需要它来处理任意矩阵，但是会有一个已知的终端状态列表和从所有其他状态到达这些状态的正概率。

目前我正在考虑使用重复乘法方法，但感觉不太理想，应该有一个函数/技巧可以在不循环的情况下计算。

我正在阅读这篇文章，但没有完全理解该方法是什么以及如何实施它。

https://math.dartmouth.edu/archive/m20x06/public_html/Lecture14.pdf

我也看过这个问题。人们似乎给出了一些手工求解的技巧，但没有给出通用算法：

https://math.stackexchange.com/questions/2003258/calculating-the-probability-of-reaching-each-absorbing-state-in-markov-chain

Answer 1

我的朋友指出了以下技巧。

特征分解意味着我们可以将原始矩阵写为

V x D x V^-1

其中 D 是具有特征值的对角矩阵，V 是特征向量。

如果我们将它自己无限次相乘，它是

V x D^inf x V^-1

我们可以使用下面的 numpy 来计算。

d, v = np.linalg.eig(m)
v @ np.diag(d >= 1).astype(int) @ np.linalg.inv(v)

由于如果对角线值