Вычислите равновесные вероятности цепи Маркова, используя итерацию Якоби в Python

Я пытаюсь вычислить функцию, которая вычисляет равновесные вероятности цепи Маркова. Для этой задачи у меня уже есть моя матрица перехода.

Теперь я пытаюсь определить функцию под названием Якоби, но не понимаю, как это сделать наиболее эффективно. Любые предложения о том, как это сделать?

До сих пор я пытался настроить его как систему уравнений и решить x = a ^ (-1) * b, но не смог правильно реализовать его из-за того, что матрица перехода является сингулярной.

Я знаю, что мне нужно умножить матрицу перехода на переменную матрицу, чтобы получить 7 отдельных уравнений. Затем мне нужно добавить уравнение x0 + x1 + x2 + x3 + x4 + x5 + x6 = 1. После того, как у меня есть все 8 уравнений, я могу решить от x0 до x6, чтобы получить свои равновесные вероятности. Вы знаете, как я могу реализовать этот процесс в коде Python?


python statistics probability markov-chains gauss
person Kalina Scarbrough    schedule 20.04.2021    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Я не уверен, что метод Якоби работает для матрицы перехода Маркова. Однако есть и другие, более простые методы нахождения стационарного распределения. Сначала решим систему уравнений, как вы описали:

import numpy as np


>>> M = np.array([[0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498], 
                  [0.5768, 0.224, 0.1494, 0.0498, 0.0, 0.0, 0.0], 
                  [0.3528, 0.224, 0.224, 0.1494, 0.0498, 0.0, 0.0], 
                  [0.1848, 0.168, 0.224, 0.224, 0.1494, 0.0498, 0.0], 
                  [0.084, 0.1008, 0.168, 0.224, 0.224, 0.1494, 0.0498]])
>>>
>>> I = np.eye(len(M)) # identity matrix
>>>
>>> A = M.T - I  # rearrange each "equation" in the system
>>>
>>> A = np.vstack([A, np.ones((1, len(A)))]) # add row of ones
>>>
>>> b = np.zeros(len(A)) # prepare solution vector
>>> b[-1] = 1.0 # add a one to the last entry
>>>
>>> np.linalg.solve(A.T @ A, A.T @ b) # solve the normal equation
array([0.22288974, 0.14776044, 0.1781388 , 0.181211  , 0.1504296 ,
       0.09080838, 0.02876204])

Вы также можете многократно умножать M само на себя, пока оно не сойдется:

>>> np.linalg.matrix_power(M, 25)[0]
array([0.22288974, 0.14776044, 0.1781388 , 0.181211  , 0.1504296 ,
       0.09080838, 0.02876204])
person blorgon    schedule 20.04.2021