R Library for Discrete Markov Chain Simulation

R library for discrete Markov chain simulation

A while back I wrote a set of functions for simulation and estimation of Discrete Markov Chain probability matrices: http://www.feferraz.net/files/lista/DTMC.R.

Relevant code for what you're asking:

simula <- function(trans,N) {
        transita <- function(char,trans) {
                sample(colnames(trans),1,prob=trans[char,])
        }

 sim <- character(N)
 sim[1] <- sample(colnames(trans),1)
 for (i in 2:N) {
  sim[i] <- transita(sim[i-1],trans)
 }

 sim
}

#example
#Obs: works for N >= 2 only. For higher order matrices just define an
#appropriate mattrans
mattrans <- matrix(c(0.97,0.03,0.01,0.99),ncol=2,byrow=TRUE)
colnames(mattrans) <- c('0','1')
row.names(mattrans) <- c('0','1')
instancia <- simula(mattrans,255) # simulates 255 steps in the process

Simulate Markov Chain using R

I am not clear what you are exactly looking for. If you are trying to obtain the transition matrix after 10 and 100 times, you can keep track of the evolution progress via Reduce, where option accumulate should be set to TRUE

P10 <- Reduce(`%*%`,replicate(10,P,simplify = FALSE),accumulate = TRUE)
P100 <- Reduce(`%*%`,replicate(100,P,simplify = FALSE),accumulate = TRUE)

markov chain simulation using R

The problem is your calculation of P1^k. In R, using ^ takes the power element by element, not through matrix multiplication. To calculate the matrix power of P1, you need to multiply it out, or use the %^% operator from the expm package. (There may be other implementations of this operation in other packages, too.)

So write your loop like this:

# Calculate probabilities for 50 steps.
P1tok <- diag(4)   # the identity matrix
for(k in 1:50){
  P1tok <- P1tok %*% P1     # multiply by another P1
  nsteps <- initState*P1tok
  one[k] <- nsteps[1,1]
  two[k] <- nsteps[1,2]
  three[k] <- nsteps[1,3]
  four[k] <- nsteps[1,4]
}

Building markov chain in r

Since you mentioned that you know how to work with statetable.msm, here's a way to translate the data into a form it can handle:

dd <- c('A-B-C-D', 'A-B-C-A', 'A-B-A-B')

Split on dashes and arrange in columns:

d2 <- data.frame(do.call(cbind,strsplit(dd,"-")))

Arrange in a data frame, identified by sequence:

d3 <- tidyr::gather(d2)

Construct the transition matrix:

statetable.msm(value,key,data=d3)

Manual simulation of Markov Chain in R (3)

I didn't inspect the rest of your code, but it seems that only prob has a mistake; you are mixing up rows with columns and instead it should be

prob <- function(simMat, fromStep, toStep, fromState, toState)
  mean(simMat[simMat[, fromStep + 1] == fromState, toStep + 1] == toState) 

Then NaN still remains a valid possibility for the following reason. We are looking at a conditional probability P(X₁=1|X₀=1) which, by definition, is well defined only when P(X₀=1)>0. The same holds with sample estimates: if there are no cases where X₀=1, then the "denominator" in the mean inside of prob is zero. Thus, it cannot and should not be fixed (i.e., returning 0 in those cases would be wrong).

R - Streamlined Markov Chain

Worked it out later: it's simplest as having a Markov chain for each age/gender group which can be simplified to a dataframe.

The initial values can be left_joined onto the transition probabilities into a data structure d.

d$temp <- lag(d$Initial * d$Terminate)
d$temp[1] <- 0 #Gets rid of NA
d$temp <- d$temp + d$hire*TotHires[1]
#where TotHires[1] represents the number hired in year 1

This gives the results after one year. For n years, we have

d$temp <- d$Initial
for (y in 1:n) {
  d$temp <- lag(d$temp * d$Terminate)
  d$temp[1] <- 0 #Gets rid of NA
  d$temp <- d$temp + d$hire*TotHires[n]
  #where TotHires[n] represents the number hired in year n
}

---

Related Topics