Search code examples
rprobabilitymarkov-chains

Simulate Steps Through a Markov Chain


I am trying to simulate a step through a Markov chain with the intent to loop through the procedure multiple times until a condition is met. (i.e., find out how many steps it takes, on average, to reach a specific state).

In this case, a state can only go one way. E.g., State 4 can transition forward to State 5, but cannot transition backward to State 3. This means the left-lower half of the transition matrix is set to zero. This is also why the method below puts arbitrarily large values in the 'prior' states. I attempt to find the correct new state by examining which probability in the specified row of the transition matrix is closest to a random number.

get_new_state <- function(current_state, trans_matrix)
{
  # generate a random number between 0-1 to compare to the transition matrix probabilities
  rand <- runif(1)
  
  # transition to where the 
  # random number falls within the transition matrix
  row <- current_state # initial condition determines the row of the trans_matrix
  col = current_state # start in the column 
  # loop thru all columns and find the correct state
  potential_states <- rep(0, each=ncol(trans_matrix)) # holds the value of the potential state it transitions to

  # in this case, we can't transition to a previous state so we set the previous state values arbitrarily high 
  # so they don't get selected in the which.min() function later
  potential_states[1:col] <- 999
  
  for(k in row:ncol(trans_matrix)) # loop thru non-zero matrix values
  {
    if(trans_matrix[row,k] > rand)
    {
      potential_states[k] <- trans_matrix[row,k] / rand
      potential_states[k] <- 1 - potential_states[k]
    }
  }
  
  # new state is equal to the index of the lowest value
  # lowest value = closest to random number
  new_state = which.min(potential_states)
  return(as.numeric(new_state))
}

I'm not sure if this approach is reasonable. I'm assuming there is a better way to simulate without the kluge that puts arbitrarily large values in potential_states[].


Solution

  • Would something like this work better (it is a one-line Markov transition):

    > new_state <- sample(1:number_states, size=1, 
    
                 prob=transition_matrix[old_state,])
    

    and just put this in a (for instance) while() loop with a counter.