概率矩阵 Numpy 的终端概率

Terminal probabilities of a probability matrix Numpy

我有一个矩阵 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.        ]]

如何在不使用重复乘法的情况下进行计算?我以为它对应于矩阵的eigenvector/eigenvalues,但是当我计算这个时我得到了一些非常不同的东西。

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 进行计算的方法?我需要它为任意矩阵工作,但会有一个已知的终端状态列表和从所有其他状态到达这些状态的正概率。

目前我正在考虑使用重复乘法,但感觉不太理想,应该有一个 function/trick 可以不用循环计算的方法。

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

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

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

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

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)

因为如果对角线值 < 1,它们将在我们相乘时趋向于 0(假设我们有一个具有有效概率的矩阵,并且所有状态都可以达到终止状态)。