Search code examples
rmarkov-chainsstochastic-process

Simulate a Markov chain path 1000 times


I have the following code for a Markov chain:

simulation.mc=function(i0,P,n.sim){
S=1:nrow(P)
X=rep(0,n.sim)
X[1]=i0
for (n in 2:n.sim){
    X[n]=sample(x=S,size=1,prob=P[X[n-1],])
    }
return(X)
}

P=matrix(
c(
    0,1/2,0,1/2,0,0,0,
    1/2,0,1/2,0,0,0,0,
    0,1,0,0,0,0,0,
    1/3,0,0,0,1/3,1/3,0,
    0,0,0,1,0,0,0,
    0,0,0,1/2,0,0,1/2,
    0,0,0,0,0,0,1
),nrow=7,byrow=T);P

X=simulation.mc(1,P,100)


T=min(which(X==7))

I have to calculate the average number of steps before reaching state 7.

I know that I need to run at least 1000 samples of the path, count the number of steps in each sample and then calculate the mean value (although some paths won't reach 7 state).

I did this, but still not working:

n.sim=100
X[i]=rep(0,n.sim)
 for (i in 1:100)
 { X[i]=simulation.mc(1,P,100)
 }

why this doesn't work? How can I include a loop inside a loop to include the function that counts the number os steps? Thanks in advance for any advice.


Solution

  • You can use replicate instead of a loop:

    replicate(1000, min(which(simulation.mc(1,P,100)==7)))
    

    @JDB provided one option for using a loop. Here are a couple more:

    # To save each entire chain in a list
    n.sim=100
    X = list()
    for (i in 1:1000) { 
      X[[i]] = simulation.mc(1,P,n.sim)
    }
    
    # To save just the number of steps to get to 7
    n.sim=100
    X = rep(NA, 1000)  
    for (i in 1:1000) { 
      X[i] = min(which(simulation.mc(1,P,n.sim)==7))
    }