Cumulative transition matrix Markov Chain simulation

Question

I have a cumulative transition matrix with probabilities for all the possible states from 1 to 5. Now the algorithm for simulating future states is following: the initial state is selected randomly, and a random value between 0 and 1 is then produced by uniform random number generator.To determine the next state in the Markov process the value of the random number is is compared with the elements of of the *i-*th row of the cumulative transition matrix determined from the preceding state. If the value of the random number is greater than the cumulative probability of the preceding state but less then or equal to the cumulative probability of the succeeding state, the succeeding state is chosen to represent the next state.

cum_trans <- matrix(c(0.1686747,0.4337349,0.6265060,0.7289157,1,
              0.2765957,0.5053191,0.6648936,0.7659574,1,
              0.2518519,0.4740741,0.6518519,0.7407407,1,
              0.1911765,0.4705882,0.6617647,0.7941176,1,
              0.2096774,0.4892473,0.6827957,0.7419355,1
              ), nrow = 5, ncol = 5, byrow = TRUE)

       [,1]      [,2]      [,3]      [,4]      [,5]
 [1,] 0.1686747 0.4337349 0.6265060 0.7289157    1
 [2,] 0.2765957 0.5053191 0.6648936 0.7659574    1
 [3,] 0.2518519 0.4740741 0.6518519 0.7407407    1
 [4,] 0.1911765 0.4705882 0.6617647 0.7941176    1
 [5,] 0.2096774 0.4892473 0.6827957 0.7419355    1


state<-matrix(NA, nrow=1,ncol=25)
i<-sample(1:5,1)
state[1]<-i
for (k in 2:25){
 rr<-runif(1)
 state[k]<-findInterval(rr,cum_trans[i,])+1
 i<-state[k]
 }
 state

So I'm not sure if I implemented the algorithm correctly. Can someone suggest modifications or improvements.

Answer 1

I see nothing wrong with your implementation, in the sense it does exactly what you described in plain English, and somewhat efficiently. Here are however a few pieces of advice, mostly about improving your coding standards.

1. Add some spaces to your code to make it more readable. Have a look at the body of basic functions (e.g. lm) to see what clean code should look like. In particular, spaces after commas and spaces on both sides of <-, =, and all binary operators.
2. Replace your hardcoded values (5 and 25) with variables. 5 can be derived from the number of rows in cum_trans and 25 is an input to your process. By assigning it to a variable (or a function input), and reusing that variable throughout your code, your code becomes easier to read and maintain, and more robust to changes. Imagine for example what would be required of you if you wanted to change the size of the transition matrix or wanted to have 50 iterations instead of 25.
3. Using a 1-row matrix is overkill and error-prone: use a vector.
4. Choose well how you name your variables and comment your code. Again, to make it easier to read and understand.
5. Learn the difference between numerics and integers. Here your states are obviously integers, yet your code stores numerics. Granted, it is not critical here, but using integers use less memory and can make your code faster. Integers are also not subject to floating point errors so they can make your code more robust in some instances.
6. Marginal speed improvement. The 25 random calls to runif(1) being independent, you can make a single call to runif(25) and store the results.
7. Think of your code in terms of inputs of outputs, then write a function.

Updated code taking some of these into account:

random.walk <- function(cum_trans, length.out = 25L) {
   num.states <- nrow(cum_trans)
   states  <- vector(mode = "integer", length = length.out)
   # pick the initial state randomly
   states[1L] <- sample(num.states, 1L)
   # jump randomly using the cum_trans matrix
   num.jumps <- length.out - 1L
   randoms <- runif(num.jumps)
   for (i in seq(num.jumps)) {
      current.state  <- states[i]
      current.prob   <- cum_trans[current.state, ]
      states[1L + i] <- findInterval(randoms[i], current.prob) + 1L
   }
   return(states)
}

Next, think of the assumptions your code is making: cum_trans must be a square matrix with each row being an increasing vector of positive values ending with 1. length.out must be an integer vector of length 1. All of these assumptions can be checked at the top of the function's body; I am only implementing a few of them here:

stopifnot(nrow(cum_trans) == ncol(cum_trans),
          all(cum_trans >= 0),
          is.integer(length.out),
          length(length.out) == 1L)

This way, there is no way you or your end-user can mess up with the inputs.

I hope you find this useful and it encourages you to write better code.