在R中绘制给定转移矩阵的马尔可夫链答案

【问题标题】：Draw markov chain given transition matrix in R在R中绘制给定转移矩阵的马尔可夫链
【发布时间】：2015-10-08 11:08:55
【问题描述】：

让trans_m 成为一阶马尔可夫链的n by n 转移矩阵。在我的问题中，n 很大，比如 10,000，而矩阵 trans_m 是由 Matrix 包构造的稀疏矩阵。否则，trans_m 的大小会很大。我的目标是在给定初始状态向量s1 和这个转移矩阵trans_m 的情况下模拟一系列马尔可夫链。考虑以下具体示例。

    n <- 5000 # there are 5,000 states in this case.
    trans_m <- Matrix(0, nr = n, nc = n, sparse = TRUE)
    K <- 5 # the maximal number of states that could be reached.
    for(i in 1:n){
        states_reachable <- sample(1:n, size = K) # randomly pick K states that can be reached with equal probability.
        trans_m[i, states_reachable] <- 1/K
    }
    s1 <- sample(1:n, size = 1000, replace = TRUE) # generate 1000 inital states
    draw_next <- function(s) {
        .s <- sample(1:n, size = 1, prob = trans_m[s, ]) # given the current state s, draw the next state .s
        .s
    }
    sapply(s1, draw_next)

给定上面的初始状态向量s1，我用sapply(s1, draw_next)来绘制下一个状态。当n 较大时，sapply 会变慢。有没有更好的办法？

【问题讨论】：

您的操作方式是获取 100 个初始状态（但您将 size=1000 放在定义 s1 的行上...）然后您只为这些状态中的每一个生成下一个状态.但是您不希望每个初始步骤都有一系列 m (m>2) 状态，而只是下一个状态，对吗？
我不明白这一点。您正在创建一个 5000*5000 矩阵，但您只使用了 5 个状态。因此，您可以将矩阵缩减为 5*5 矩阵以进行模拟。
@ColonelBeauvel 我原来的帖子有错字。我打算画 1000 个初始状态，尽管我在评论中输入了 100。感谢您指出。我无法将转换矩阵减少到 5 乘 5，因为五个可达到的下一个状态相对于条件变量（当前状态）会发生变化。

标签： r matrix markov-chains

【解决方案1】：

按行重复索引可能会很慢，因此处理转换矩阵的转置并使用列索引并从内部函数中分解索引会更快：

R>    trans_m_t <- t(trans_m)
R>
R>    require(microbenchmark)
R>    microbenchmark(
+       apply(trans_m_t[,s1], 2,sample, x=n, size=1, replace=F)
+     ,
+       sapply(s1, draw_next)
+     )
Unit: milliseconds
                                                            expr        min
 apply(trans_m_t[, s1], 2, sample, x = n, size = 1, replace = F) 111.828814
                                           sapply(s1, draw_next) 499.255402
          lq        mean      median          uq        max neval
 193.1139810 190.4379185 194.6563380 196.4273105 270.418189   100
 503.7398805 512.0849013 506.9467125 516.6082480 586.762573   100

由于您已经在使用稀疏矩阵，因此您也许可以通过直接处理三元组获得更好的性能。使用更高级别的矩阵运算符可以触发重新压缩。

【讨论】：

转置是个好主意。不过我试过三胞胎。这非常方便。谢谢。