转移概率的递归函数（Chapman-Kolmogorov eq.）答案

【问题标题】：Recursive function (Chapman-Kolmogorov eq.) for Transition Probabilities转移概率的递归函数（Chapman-Kolmogorov eq.）
【发布时间】：2018-03-19 11:55:27
【问题描述】：

我一直在构建一个递归函数，最好通过一个简单的例子来说明。

采用具有 2 个状态，状态 1 和状态 2 的马尔科夫过程。符号 p_ij 表示在当前状态为 的情况下转换到状态 j 的概率>我。在这个例子中，

p_11 = 0.8（在当前状态为状态 1 的情况下停留在状态 1 的概率）
p_12 = 0.2
p_21 = 0.6
p_22 = 0.4

而转移概率矩阵为：

import numpy as np
pij = np.array([[.8, .2], [.6, .4]])
print(pij)
# [[ 0.8  0.2]
#  [ 0.6  0.4]]

n 步转移概率，表示为 r_ij(n)，表示在 n 个时间段之后的状态为 n 的概率em>j，假设当前状态是i。 r_ij(n) 可以使用 Chapman-Kolmogorov 方程找到，

有初始条件

m 是状态的总数。 [来自 Bertsekas/Tsitsiklis，2008 年。]

我正在尝试构建 r_ij(n)。前 5 个步骤应如下所示：

我的开始：

def r(p, n):
    m = np.sqrt(p.size)  # or p.shape[0]
    if n == 1:
        return p
    elif n > 1:
        res = []
        for k in range(m):
            for i in p:
                for j in i:
                    # This line is patently wrong...
                    # Not sure how to reference i
                    return r(n - 1) * p[k, j]
        return np.sum(res)

p0 = np.array([[.8, .2], [.6, .4]])
print(r(p0, n=5))
# [[.7501, .2499],
#  [.7498, .2502]]

但是我对这个符号有点迷茫。

【问题讨论】：

你已经有了一个好的开始。当您定义pij = np.array([[.8, .2], [.6, .4]]) 时，更自然的标识符应该是p。类似地，rij 应该是 r。您的定义说 3 个循环：在 i、j、k 上，但您只在 k 上编写了一个循环。请注意，您使用的是单下标，其中双下标（如p[k, j]）是合适的。

标签： python numpy recursion

【解决方案1】：

您不一定需要递归函数。 matrix_power 本质上就是你要找的rij：

def rij(pij, n):
    return matrix_power(pij, n)

pij = np.array([[.8, .2], [.6, .4]])

from numpy.linalg.linalg import matrix_power

matrix_power(pij, 2)
#array([[ 0.76,  0.24],
#       [ 0.72,  0.28]])

matrix_power(pij, 3)
#array([[ 0.752,  0.248],
#       [ 0.744,  0.256]])

matrix_power(pij, 4)
#array([[ 0.7504,  0.2496],
#       [ 0.7488,  0.2512]])

matrix_power(pij, 5)
#array([[ 0.75008,  0.24992],
#       [ 0.74976,  0.25024]])

要定义递归函数，np.dot 将使任务更容易：

def rij(pij, n):
    if n == 1:
        return pij
    else:
        return np.dot(rij(pij, n-1), pij)

rij(pij, 5)
#array([[ 0.75008,  0.24992],
#       [ 0.74976,  0.25024]])

【讨论】：

很棒的答案。使用dot，因为矩阵幂实际上是重复矩阵平方+乘法？
是的。这基本上就是matrix_power 所做的，重复矩阵乘法。使用带有递归函数的np.dot来模拟这个过程。
对于n 是 NumPy 数组的情况，我有什么办法可以矢量化？ n == 1 测试不喜欢这种情况，matrix_power 需要 n 作为 int。结果将是 3d，形状 (len(n), m, m)。
给@Psidom 的答案加点，matrix_power 比直接递归乘法更快，因为该函数使用动态编程来减少所需的乘法次数。
我认为不会有有效的矢量化。如果您要为多个ns 计算此值，我建议在此处介绍一些记忆技术。例如，对于[10, 30, 40] 的幂，您可能不想从地面计算所有内容，而是想根据10 计算30，例如基于40 的30。