График зависимости вероятности P (Xn = s) от времени?

Question

Вся проблема: у нас есть модель цепей Маркова с 5 состояниями: s, t, m, f, r

TPM:

      P <- matrix(c(.84,.03,.01,.03,.03,
                    .11,.80,.15,.19,.09,
                    .01,.04,.70,.02,.05,
                    .04,.10,.07,.75,.00,
                    .00,.03,.07,.01,.83),
                  nrow=5

                   )

С матричным умножением, ограничение Распределение выходит:

(0.1478365,0.4149259,0.09555939,0.2163813,0.1252968)

Я пытаюсь построить график изменения P (Xn = s) как функции времени.

Учитывая, что начальное распределение P (X0 = i) = 1/5, т.е.:

                  s   t   m   f   r  
           α = ( 1/5 1/5 1/5 1/5 1/5 )

Мне нужно построить график P (Xn = s) (по оси y) относительно n = 0 , 1, 2, 3, 4, 5 (ось X).

Answer 1

Вот базовое решение R, которое позволяет работать на матричной мощности и просматривать ход эволюции

n <- 20
res <- do.call(rbind,Reduce(`%*%`,c(list(a),replicate(n,P,simplify = FALSE)),accumulate = T))

затем, аналогично подходу построения графиков @ 42- , вы можно использовать, чтобы увидеть эволюцию относительно 1:n

matplot(seq(nrow(r)),r)

Answer 2

Пакет expm имеет функцию мощности матрицы, %^%. Сначала я проверил 10-е состояние, но оно не достигло стабильности, поэтому я пошел с 20-м состоянием

Каждое последующее состояние находится в прогрессии a %*% P, a %*% (P%^%2), ...., a %*% (P%^%20)

library(expm)

#Attaching package: ‘expm’

#The following object is masked from ‘package:Matrix’:

 #   expm

Так что, думаю, у меня было это уже доступно, неважно.

evolve <- sapply(1:20, function(n){ a%*% (P %^% n)})  # 5 x 20 matrix
png(); matplot( 1:20, t(evolve) );dev.off()  # matplot needs data in rows