从给定序列构建N阶马尔可夫转移矩阵

Question

从给定序列构建N阶马尔可夫转移矩阵

3

我正在尝试创建一个函数，可以将给定的输入序列转换为所请求的转移矩阵。我找到了第一阶马尔可夫转移矩阵的实现。

现在，我希望能够想出一个解决方案，可以计算出第二和第三阶转移矩阵。

以下是第一阶矩阵实现的示例：

import numpy as np

# sequence with 3 states -> 0, 1, 2

a = [0, 1, 0, 0, 0, 2, 2, 1, 1, 1, 0, 0, 0, 0, 0, 1, 2, 2, 2, 0, 0, 2]


def transition_matrix_first_order(seq):
    M = np.full((3, 3), fill_value = 1/3, dtype= np.float64)
    for (i,j) in zip(seq, seq[1:]):
        M[i, j] += 1

    M = M / M.sum(axis = 1, keepdims = True)

    return M

print(transition_matrix_first_order(a))

这给了我这个：

[[0.61111111 0.19444444 0.19444444]
 [0.38888889 0.38888889 0.22222222]
 [0.22222222 0.22222222 0.55555556]]

当制作二阶矩阵时，它应该有unique_state_count ** order行和unique_state_count列。在上面的例子中，我有3个唯一状态，因此矩阵将具有9x3的结构。理想的函数示例： cal_tr_matrix(seq, unique_state_count, order)

- Emil Mirzayev

当你增加矩阵时，你应该简单地考虑到最后一个先前的索引。此外，我会将M = M/M.sum(..)从for循环中移除，因为这看起来不正确。 - undefined

抱歉，我会进行编辑的。我在从笔记本上复制粘贴时犯了一个错误 :) - undefined

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Richard Nemeth · Accepted Answer

我认为你对马尔可夫链及其转移矩阵有一些误解。

首先，你的函数生成的估计转移矩阵是不正确的。为什么呢？我们来复习一下。

离散时间内有N个不同状态的离散马尔可夫链具有大小为N x N的转移矩阵P，其中（i，j）元素是P（X_1 = j | X_0 = i），即从状态i在单个时间步骤中转移到状态j的概率。

现在，阶数为n的转移矩阵，表示为P^{n}，再次是大小为N x N的矩阵，其中（i，j）元素是P（X_n = j | X_0 = i），即在n个时间步骤中从状态i转移到状态j的概率。

一个美妙的结果表明：P^{n} = P^n，即将单步转移矩阵的n个幂次相加会得到n步转移矩阵。

现在，通过这个复习，需要做的就是从给定序列中估计P，然后可以使用已经估计出的P来估算P^{n}，并将矩阵取n次幂。那么如何估计矩阵P呢？如果我们将N_{ij}表示为从状态i到状态j的转移观察次数，将N_{i*}表示为处于状态i的观察次数，则P_{ij} = N_{ij} / N_{i*}。

最终Python代码如下：

import numpy as np

def transition_matrix(arr, n=1):
    """"
    Computes the transition matrix from Markov chain sequence of order `n`.

    :param arr: Discrete Markov chain state sequence in discrete time with states in 0, ..., N
    :param n: Transition order
    """

    M = np.zeros(shape=(max(arr) + 1, max(arr) + 1))
    for (i, j) in zip(arr, arr[1:]):
        M[i, j] += 1

    T = (M.T / M.sum(axis=1)).T

    return np.linalg.matrix_power(T, n)

transition_matrix(arr=a, n=1)

>>> array([[0.63636364, 0.18181818, 0.18181818],
>>>       [0.4       , 0.4       , 0.2       ],
>>>       [0.2       , 0.2       , 0.6       ]])

transition_matrix(arr=a, n=2)

>>> array([[0.51404959, 0.22479339, 0.26115702],
>>>       [0.45454545, 0.27272727, 0.27272727],
>>>       [0.32727273, 0.23636364, 0.43636364]])

transition_matrix(arr=a, n=3)

>>> array([[0.46927122, 0.23561232, 0.29511645],
>>>       [0.45289256, 0.24628099, 0.30082645],
>>>       [0.39008264, 0.24132231, 0.36859504]])

有趣的是，当你将顺序n设置为一个相当高的数字时，矩阵P的高次幂似乎会收敛到一些非常特定的值。这被称为马尔可夫链的稳态/不变分布，并且它很好地说明了链在长时间/转换期内的行为。此外：

P = transition_matrix(a, 1)
P111 = transition_matrix(a, 111)
print(P)
print(P111.dot(P))

编辑：根据您的评论，我建议使用更高维度的矩阵来表示更高阶的情况，而不是增加行数。一种方法是这样的：

def cal_tr_matrix(arr, order):

    _shape = (max(arr) + 1,) * (order + 1)
    M = np.zeros(_shape)

    for _ind in zip(*[arr[_x:] for _x in range(order + 1)]):
        M[_ind] += 1
    return M

res1 = cal_tr_matrix(a, 1)
res2 = cal_tr_matrix(a, 2)

现在，res1[i, j]元素表示发生了多少次i->j的转换，而res2[i, j, k]元素表示发生了多少次i->j->k的转换。