моделирование цепей Маркова с использованием R

Question

Я хочу построить график распределения вероятностей состояний в разное время t.У меня есть матрица перехода

P1 <- matrix(c(0, 1, 0, 0, 0, 0, 2/3, 1/3, 0, 1, 0, 0, 0, 0, 0, 1), 4, 4, byrow=TRUE) 

Чтобы построить график от p1 до p50, я использую следующий код.

library(ggplot2)

initState <- c(1,0,0,0) # The initial state will be state 1

# Initiate probability vectors
one <- c()
two <- c()
three <- c()
four <- c()

# Calculate probabilities for 50 steps.
for(k in 1:50){
  nsteps <- initState*P1^k
  one[k] <- nsteps[1,1]
  two[k] <- nsteps[1,2]
  three[k] <- nsteps[1,3]
  four[k] <- nsteps[1,4]
}

# Make dataframes and merge them
one <- as.data.frame(one)
one$Group <- 'one'
one$Iter <- 1:50
names(one)[1] <- 'Value'

two <- as.data.frame(two)
two$Group <- 'two'
two$Iter <- 1:50
names(two)[1] <- 'Value'

three <- as.data.frame(three)
three$Group <- 'three'
three$Iter <- 1:50
names(three)[1] <- 'Value'

four <- as.data.frame(four)
four$Group <- 'four'
four$Iter <- 1:50
names(four)[1] <- 'Value'

steps <- rbind(one,two,three,four)

# Plot the probabilities using ggplot
ggplot(steps, aes(x = Iter, y = Value, col = Group))+
  geom_line() +
  xlab('Chain Step') +
  ylab('Probability') +
  ggtitle('50 Step Chain Probability Prediction')+
  theme(plot.title = element_text(hjust = 0.5))

Но результат очень странный, кто-нибудь может узнать, где я ошибся в своем коде?

Answer 1

Проблема в том, что вы вычислили P1^k.В R использование ^ берет степенной элемент за элементом, а не через матричное умножение.Чтобы рассчитать мощность матрицы P1, вам нужно умножить ее или использовать оператор %^% из пакета expm.(Могут быть и другие реализации этой операции в других пакетах.)

Итак, напишите свой цикл так:

# Calculate probabilities for 50 steps.
P1tok <- diag(4)   # the identity matrix
for(k in 1:50){
  P1tok <- P1tok %*% P1     # multiply by another P1
  nsteps <- initState*P1tok
  one[k] <- nsteps[1,1]
  two[k] <- nsteps[1,2]
  three[k] <- nsteps[1,3]
  four[k] <- nsteps[1,4]
}