Построение матрицы марковских переходов N-го порядка из заданной последовательности - PullRequest
0 голосов
/ 22 сентября 2019

Я пытаюсь создать функцию, которая может преобразовать заданную входную последовательность в матрицу перехода запрошенного порядка.Я нашел реализацию для матрицы марковских переходов первого порядка.

Теперь я хочу найти решение, которое может вычислять матрицы переходов 2-го и 3-го порядка.

Примерреализация матрицы 1-го порядка:

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]]

При создании матрицы 2-го порядка она должна иметь unique_state_count ** order строк и unique_state_count столбцов.В приведенном выше примере у меня есть 3 уникальных состояния, поэтому матрица будет иметь структуру 9x3.Желаемый образец функции: cal_tr_matrix(seq, unique_state_count, order)

1 Ответ

0 голосов
/ 22 сентября 2019

Я думаю, что у вас есть небольшое недопонимание о цепях Маркова и их матрицах переходов.

Прежде всего, оценочная матрица переходов, которую создает ваша функция, к сожалению, не верна.Почему?Давайте освежимся.

Дискретная цепь Маркова в дискретном времени с 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->К случилось.

...