如何在 R 中制作第二个马尔可夫链模型或更高模型的可重现示例？

Question

我将开始定义一个值向量，如下所示：

data <- c(2.3, 1.7, 1.7, 1.4, 1.7, 2.0, 2.2, 1.8, 1.7, 1.6, 2.0, 1.5, 1.1, 1.4, 1.8, 2.0, 1.5, 1.2, 1.0, 1.2, 1.5, 1.6, 1.1,
      1.5, 2.0, 2.1, 2.1, 2.0, 1.7, 1.7, 1.7, 1.3, 0.8, -0.1, 0.0, -0.1, -0.2, 0.0, 0.1, 0.2, 0.2, 0.0, 0.2, 0.5, 0.7, 1.4,
      1.0, 0.9, 1.1, 1.0, 1.0, 0.8, 1.1, 1.5, 1.6, 1.7, 2.1, 2.5, 2.7, 2.4, 2.2, 1.9, 1.6, 1.7, 1.9, 2.2, 2.0, 2.2, 2.1,
      2.1, 2.2, 2.4, 2.5, 2.8, 2.9, 2.9, 2.7, 2.3, 2.5, 2.2, 1.9, 1.6, 1.5, 1.9, 2.0, 1.8, 1.6, 1.8, 1.7, 1.7, 1.8, 2.1,
      2.3, 2.5, 2.3, 1.5, 0.3, 0.1, 0.6, 1.0, 1.3, 1.4, 1.2, 1.2, 1.4, 1.4, 1.7, 2.6, 4.2, 5.0, 5.4, 5.4, 5.3, 5.4, 6.2,
      6.8, 7.0, 7.5, 7.9, 8.5, 8.3, 8.6, 9.1, 8.5, 8.3, 8.2, 7.7, 7.1, 6.5, 6.4, 6.0, 5.0, 4.9, 4.0, 3.0, 3.2, 3.7, 3.7)

从这个向量中，我做了一个 ifelse 语句来表示它的“增加”、“减少”和“停滞”

direction <- c(0,diff(data))
direction_states <- ifelse(direction > 0, "Increase", ifelse(direction < 0, "Decrease", "Stagnant"))

最后从这两个相应的键中得到一个数据框：

df_track <- as.data.frame(cbind(data, direction_states))
df_track <- df_track[-1,]

这就是我使用包构建一阶马尔可夫链的方式

markovchain

library(markovchain)

transition_matrix <- as.matrix(table(df_track$direction_states, lag(df_track$direction_states)) / rowSums(table(df_track$direction_states, lag(df_track$direction_states))))
transition_matrix <- as.numeric(transition_matrix)
dim(transition_matrix) <- c(3, 3)
transition_matrix[is.na(transition_matrix)] <- 0
markov_chain <- new("markovchain", states = c("Decrease", "Increase", "Stagnant"), 
                    transitionMatrix = transition_matrix, name = "Markov 1st Order Value")
print(markov_chain)
set.seed(123) 
simulated_states <- rmarkovchain(n = 10, object = markov_chain, t0 = tail(df_track$direction_states, 1))
print(simulated_states)

但是我无法重现制作第二个或更高马尔可夫模型的示例，因为以下陈述：

markov_chain <- new("markovchain", states = c("Decrease", "Increase", "Stagnant"), 
                        transitionMatrix = transition_matrix, name = "Markov 1st Order Value")

只允许输入 2D 矩阵，不允许输入 3D 或以上。显示这样的错误：

Error in validObject(.Object) : 
  invalid class “markovchain” object: invalid object for slot "transitionMatrix" in class "markovchain": got class "array", should be or extend class "matrix"

感谢帮助！

Answer 1

我想我刚刚找到了答案，当我尝试像这样执行它时，在

fitHigherOrder()

包中有一个名为

markovchain

的函数：

high_order4 <- fitHigherOrder(df_track$direction_states, 4) #using 4 orders

我得到这些结果：

$lambda
[1] 2.026152e-08 3.087294e-07 9.999995e-01 1.391219e-07

$Q
$Q[[1]]
          Decrease   Increase   Stagnant
Decrease 0.5614035 0.25373134 0.58333333
Increase 0.3333333 0.65671642 0.33333333
Stagnant 0.1052632 0.08955224 0.08333333

$Q[[2]]
           Decrease  Increase   Stagnant
Decrease 0.47368421 0.3787879 0.33333333
Increase 0.49122807 0.4848485 0.58333333
Stagnant 0.03508772 0.1363636 0.08333333

$Q[[3]]
          Decrease   Increase   Stagnant
Decrease 0.4035088 0.46153846 0.16666667
Increase 0.4912281 0.46153846 0.75000000
Stagnant 0.1052632 0.07692308 0.08333333

$Q[[4]]
           Decrease   Increase   Stagnant
Decrease 0.42857143 0.44615385 0.16666667
Increase 0.48214286 0.46153846 0.75000000
Stagnant 0.08928571 0.09230769 0.08333333


$X
  Decrease   Increase   Stagnant 
0.41605839 0.48905109 0.09489051 

Q[[1]]

与上面的一阶马尔可夫链的结果几乎相同，除了这个打印了列式结果（总和为 1 的概率），而上面的一阶马尔可夫链打印了行式结果。因此，如果我这样做，我可以重现它：

for(g in 1:length(order4$Q)){
  order4$Q[[g]] <- t(order4$Q[[g]])
}

markov_chain_order2 <- new("markovchain", states = c("Decrease", "Increase", "Stagnant"), 
                    transitionMatrix = order4$Q[[2]], name = "Markov 2nd Order Value")
print(markov_chain_order2)
set.seed(123) # for reproducibility
simulated_states <- rmarkovchain(n = 10, object = markov_chain_order2, t0 = tail(df_track$direction_states, 1))
print(simulated_states)

如果 Q[[2]] 及以上表示其在 2D 矩阵中的二阶马尔可夫链概率或以上。 （不过，如果我错了，请纠正我）