离散马尔可夫链模拟的R语言库

Question

我正在寻找类似于“msm”包的东西，但是针对离散马尔可夫链。例如，如果我有一个定义为下面这样的转移矩阵：

Pi <- matrix(c(1/3,1/3,1/3,
0,2/3,1/6,
2/3,0,1/2))

针对状态A、B、C，如何根据该转移矩阵模拟马尔可夫链？

- stevejb

---

NVM：找到答案了http://books.google.com/books?id=AALhk_mt7SYC&lpg=PA116&dq=r%20markov%20chain&pg=PA119#v=onepage&q=r%20markov%20chain&f=false - stevejb

Answer 1

不久前，我编写了一组用于模拟和估计离散马尔可夫链概率矩阵的函数：http://www.feferraz.net/files/lista/DTMC.R。

以下是您所需的相关代码：

simula <- function(trans,N) {
        transita <- function(char,trans) {
                sample(colnames(trans),1,prob=trans[char,])
        }

 sim <- character(N)
 sim[1] <- sample(colnames(trans),1)
 for (i in 2:N) {
  sim[i] <- transita(sim[i-1],trans)
 }

 sim
}

#example
#Obs: works for N >= 2 only. For higher order matrices just define an
#appropriate mattrans
mattrans <- matrix(c(0.97,0.03,0.01,0.99),ncol=2,byrow=TRUE)
colnames(mattrans) <- c('0','1')
row.names(mattrans) <- c('0','1')
instancia <- simula(mattrans,255) # simulates 255 steps in the process

- Fernando H Rosa

Answer 2

啊，我正在为你写解释的时候，你已经找到了答案。这里是我想出来的一个简单示例：

run = function()
{
    # The probability transition matrix
    trans = matrix(c(1/3,1/3,1/3,
                0,2/3,1/3,
                2/3,0,1/3), ncol=3, byrow=TRUE);

    # The state that we're starting in
    state = ceiling(3 * runif(1, 0, 1));
    cat("Starting state:", state, "\n");

    # Make twenty steps through the markov chain
    for (i in 1:20)
    {
        p = 0;
        u = runif(1, 0, 1);

        cat("> Dist:", paste(round(c(trans[state,]), 2)), "\n");
        cat("> Prob:", u, "\n");

        newState = state;
        for (j in 1:ncol(trans))
        {
            p = p + trans[state, j];
            if (p >= u)
            {
                newState = j;
                break;
            }
        }

        cat("*", state, "->", newState, "\n");
        state = newState;
    }
}

run();

请注意，您的概率转移矩阵每行的总和不为1，它应该是1。我举的例子使用了稍微修改过的概率转移矩阵，符合这个规则。

- icio

---

感谢您的回答。您的代码非常易读，我真的很感激。 - stevejb

可读性强的代码？在我的经验中，这个概念已经完全失去了大多数使用 R 的人的注意。希望这有所帮助！ - icio

为了生成1至3之间的随机整数，我认为使用sample(1:3, 1)比使用ceiling(3 * runif(1, 0, 1))更容易理解。对于最内层的for循环，你可以简单地使用newState <- sample(1:ncol(trans), 1, prob=trans[state,])。这更清晰地展示了正在发生的事情。而且，你甚至不必再归一化行了。 - Ken Williams

Answer 3

您现在可以使用CRAN中可用的markovchain包。 用户手册非常好，有几个示例。

- Rodrigo Zepeda