在R中给定转移矩阵，绘制马尔可夫链。

Question

在R中给定转移矩阵，绘制马尔可夫链。

4

让trans_m成为一个n乘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变得很慢。有更好的方法吗？

- semibruin

你的操作方式是，你取100个初始状态（但在定义s1的那行代码中加上了size=1000...），然后只为每个状态生成下一个状态。但你不想为每个初始步骤生成一系列m（m>2）个状态，而只需要下一个状态，对吗？ - Colonel Beauvel

我不明白你的意思。你创建了一个5000 * 5000的矩阵，但是你只使用了5个状态。因此，你可以将矩阵缩小为5 * 5的矩阵用于你的模拟。 - Colonel Beauvel

@ColonelBeauvel 在我原始帖子中有一个拼写错误。我本打算绘制1000个初始状态，但在评论中却输入了100。谢谢你指出来。我无法将转移矩阵缩小到5乘5的大小，因为五个可达到的下一个状态会随着条件变量（当前状态）而改变。 - semibruin

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Neal Fultz · Accepted Answer

重复按行索引可能会很慢，因此更快的方法是使用转置的转移矩阵进行列索引，并将索引从内部函数中分解出来：

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

既然您已经在使用稀疏矩阵，那么直接在三元组上操作可能会获得更好的性能。使用更高级的矩阵运算符可能会触发重新压缩。