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[]
.
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.