Я думаю, что у вас есть небольшое недопонимание о цепях Маркова и их матрицах переходов.
Прежде всего, оценочная матрица переходов, которую создает ваша функция, к сожалению, не верна.Почему?Давайте освежимся.
Дискретная цепь Маркова в дискретном времени с N
различными состояниями имеет матрицу перехода P
размера N x N
, где элемент (i, j) равен P(X_1=j|X_0=i)
, то есть вероятность переходаиз состояния i
в состояние j
за один шаг по времени.
Теперь матрица перехода порядка n
, обозначенная P^{n}
, снова представляет собой матрицу размера N x N
, где a (i, j) элемент равен P(X_n=j|X_0=i)
, т. е. вероятность перехода из состояния i
в состояние j
за n
шагов по времени.
Прекрасный результат говорит: P^{n} = P^n
, т.е. принятие n
Мощность матрицы одношагового перехода дает вам матрицу n
-шагового перехода.
Теперь с этим резюме все, что нужно, это оценить P
из заданной последовательности, а затем оценить P^{n}
можно просто использовать уже оцененный P
и взять 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->К случилось.