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.
- 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 aroundafter commas and spaces on both sides of <-
, =
, and all binary operators.
- 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.
- Using a 1-row matrix is overkill and error-prone: use a vector.
- Choose well how you name your variables and comment your code. Again, to make it easier to read and understand.
- 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.
- 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.
- Think of your code in terms of inputs of outputs, then write a function.
Updated code taking some of these into account:
numrandom.stateswalk <- nrowfunction(cum_trans, length.out = 25L) {
num.jumps states <- 24Lnrow(cum_trans)
states <- vector(modelmode = "integer", length = 1L + numlength.jumpsout)
randoms <- runif(num.jumps)
# 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_jumpsnum.jumps)) {
current.state <- states[i]
current.prob <- cum_trans[current.state, ]
states[1L + i] <- findInterval(randoms[i], current.prob) + 1L
}
return(states)
}