3

I have a sparse transition matrix (more specifically an object of class dgCMatrix, which is a class defined in Matrix package for which rowSums is always equal to 1 ) and I'd like to use it as a transition matrix to generate a Markov Chain.

Does anyone know about a function (in any package, or user built) to sample a discrete Markov Chain using that object as transition matrix?

I've already looked at the markovchain package, but for the functions in that package it seems that is required a traditional (dense) matrix.

Thanks

edit: here it is my quick and dirty solution:

simulate <- function(TransMat,N) {

  # idx is in 1:length(TransMat@Dimnames[[1]]), 
  # while the values in TransMat@i are in 0:(length(TransMat@Dimnames[[1]]) - 1)

  transition <- function(idx,TransMat) {
    # Returns the index (in 1:length(TransMat@Dimnames[[1]])) of the next state,
    # randomly sampled according to the idx-th row of the transition matrix TransMat.

    # Discrete variable sampling
    u = runif(1)                                       
    cs = cumsum(TransMat@x[(TransMat@i + 1) == idx])   
    # Position of the sampled value inside the vector of indices with value equal to idx, TransMat@i[(TransMat@i + 1) == idx]
    # which(cs > u)[1]
    # Position of the sampled value inside whole the vector of indices, TransMat@i
    # which((TransMat@i + 1) == idx)[which(cs > u)[1]]
    which(TransMat@p >= which((TransMat@i + 1) == idx)[which(cs > u)[1]])[1] - 1
  }

  sim <- integer(N)
  # The initial state is fixed
  sim[1] <- 1
  # The initial state is randomly distributed according to a given (discrete) initial distribution:
  # TO BE DONE
  for (i in 2:N) 
    sim[i] <- transition(sim[i-1], TransMat)

  sim
}
NDonelli
  • 51
  • 6
  • this shouldn't be too hard, via `?rmultinom`, but you might want to rephrase the question as "how can I program this?" rather than "is there some existing code to do this?" – Ben Bolker Sep 20 '16 at 18:39
  • Not exactly. I know how to code it (although I admit I haven't already done it!). I'd like to know if there is an already implemented solution which almost surely would be more efficient wrt the solution I could implement quickly. Anyway, thanks a lot for the answer! – NDonelli Sep 20 '16 at 18:47
  • maybe post some code and ask for improvements? That would give people a way to suggest alternatives (although I am not optimistic that one exists) *or* ways of making your code efficient (although such a question might also be recommended for migration to Code Review SE ...) – Ben Bolker Sep 20 '16 at 18:52
  • OK, I am gonna do it asap (maybe directly on Code Review SE...). Thanks a lot! – NDonelli Sep 20 '16 at 19:00
  • 2
    As promised, I just added my quick and dirty solution (loosely based on a solution for a similar problem proposed in this question [link](http://stackoverflow.com/questions/2754469/r-library-for-discrete-markov-chain-simulation) by Fernando H Rosa. Any comment or improvement? – NDonelli Sep 21 '16 at 09:20

0 Answers0