7

I am looking for something like the 'msm' package, but for discrete Markov chains. For example, if I had a transition matrix defined as such

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

for states A,B,C. How can I simulate a Markov chain according to that transition matrix?

1

3 Answers 3

8

A while back I wrote a set of functions for simulation and estimation of Discrete Markov Chain probability matrices: http://www.feferraz.net/files/lista/DTMC.R.

Relevant code for what you're asking:

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
6

Argh, you found the solution whilst I was writing it up for you. Here's a simple example that I came up with:

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();

Note that your probability transition matrix doesn't sum to 1 in each row, which it should do. My example has a slightly altered probability transition matrix which adheres to this rule.

3
  • Thanks for the answer. Your code is very readable. I really appreciate it.
    – stevejb
    May 2, 2010 at 20:56
  • 2
    The readable code? In my experience this concept has been totally lost on most of the people who use R. Hope it helps!
    – icio
    May 2, 2010 at 22:09
  • 4
    To generate a random integer from 1 to 3, I think sample(1:3, 1) is a bit easier to grok than ceiling(3 * runif(1, 0, 1)). Also, for the innermost for-loop, you can simply use newState <- sample(1:ncol(trans), 1, prob=trans[state,]). That shows more clearly what's going on. And then you won't even have to normalize the rows, either. May 3, 2010 at 14:51
6

You can now use the markovchain package available in CRAN. The user manual. is pretty good and has several examples.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.