【问题标题】：Python code explanation for stationary distribution of a Markov chain马尔可夫链平稳分布的Python代码解释
【发布时间】：2022-05-09 09:56:10
【问题描述】：

我有这个代码：

import numpy as np
from scipy.linalg import eig 
transition_mat = np.matrix([
    [.95, .05, 0., 0.],\
    [0., 0.9, 0.09, 0.01],\
    [0., 0.05, 0.9, 0.05],\
    [0.8, 0., 0.05, 0.15]])

S, U = eig(transition_mat.T)
stationary = np.array(U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat)
stationary = stationary / np.sum(stationary)

>>> print stationary
[ 0.34782609  0.32608696  0.30434783  0.02173913]

但我看不懂这行：

stationary = np.array(U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat)

谁能解释一下：U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat？

我知道例程返回S: eigenvalue, U: eigenvector。我需要找到特征值1对应的特征向量。我写了下面的代码：

for i in range(len(S)):
    if S[i] == 1.0:
        j = i

 matrix = np.array(U[:, j].flat)

我得到了输出：

:  [ 0.6144763   0.57607153  0.53766676  0.03840477]

但它没有给出相同的输出。为什么？！

【问题讨论】：

嗯...它看起来从 U 中抓取元素，其中 S 中的相应元素的值为 1...但我不确定那是什么[0][0] 是为了。然后它将整个事物展平为一维数组。
嗨，凯文，感谢您的回答。你能告诉我这段代码是如何工作的：where(np.abs(S-1.)
你试过一个简单的例子吗？阅读任何文档？ 究竟是什么让你感到困惑？
@jonrsharpe，是的，我还在尝试。我已经用我的代码更新了这个问题。请检查。谢谢。
这里的页边距太小，无法完整地写下所有内容。只需仔细阅读numpy.where 的文档，然后尝试示例！ ;)

标签： python python-3.x numpy

【解决方案1】：

如何找到平稳分布。

好的，我来这篇文章是想看看是否有内置的方法可以找到平稳分布。好像没有。因此，对于来自 Google 的任何人，这就是我在这种情况下找到平稳分布的方法：

import numpy as np

#note: the matrix is row stochastic.
#A markov chain transition will correspond to left multiplying by a row vector.
Q = np.array([
    [.95, .05, 0., 0.],
    [0., 0.9, 0.09, 0.01],
    [0., 0.05, 0.9, 0.05],
    [0.8, 0., 0.05, 0.15]])

#We have to transpose so that Markov transitions correspond to right multiplying by a column vector.  np.linalg.eig finds right eigenvectors.
evals, evecs = np.linalg.eig(Q.T)
evec1 = evecs[:,np.isclose(evals, 1)]

#Since np.isclose will return an array, we've indexed with an array
#so we still have our 2nd axis.  Get rid of it, since it's only size 1.
evec1 = evec1[:,0]

stationary = evec1 / evec1.sum()

#eigs finds complex eigenvalues and eigenvectors, so you'll want the real part.
stationary = stationary.real

那条奇怪的线在做什么。

让我们把这条线分成几部分：

#Find the eigenvalues that are really close to 1.
eval_close_to_1 = np.abs(S-1.) < 1e-8

#Find the indices of the eigenvalues that are close to 1.
indices = np.where(eval_close_to_1)

#np.where acts weirdly.  In this case it returns a 1-tuple with an array of size 1 in it.
the_array = indices[0]
index = the_array[0]

#Now we have the index of the eigenvector with eigenvalue 1.
stationary = U[:, index]

#For some really weird reason, the person that wrote the code
#also does this step, which is completely redundant.
#It just flattens the array, but the array is already 1-d.
stationary = np.array(stationary.flat)

如果将所有这些代码行压缩成一行，则会得到stationary = np.array(U[:, np.where(np.abs(S-1.)<1e-8)[0][0]].flat)

如果你删除多余的东西，你会得到stationary = U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]]

为什么您的代码会给出不同的固定向量。

正如@Forzaa 指出的那样，您的向量不能表示概率向量，因为它的总和不为 1。如果将它除以它的总和，您将得到原始代码 sn-p 具有的向量。

只需添加这一行：

stationary = matrix/matrix.sum()

然后您的固定分布将匹配。

【讨论】：

【解决方案2】：

stationary = np.array(U[:,np.where(np.abs(S-1.) < 1e-8)[0][0]].flat)

这段代码在U中寻找对应特征值的元素——1小于1e-8

【讨论】：

【解决方案3】：

实际上，只需做一个简单的 while 迭代。我以随机 P 为例

def get_stationary(n):
    row = n
    pi = np.full((1, row), 1 / row)
    T = np.array([[1/4,1/2,1/4],
                  [1/3,0,2/3],
                  [1/2,0,1/2]])
    while True:
        new_pi = np.dot(pi, T)
        if np.allclose(pi, new_pi):
            return pi 
            break
        pi = new_pi
print(get_stationary(3))

【讨论】：