1

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.

Maruska
  • 13
  • 4

2 Answers2

1

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))
}
eipi10
  • 91,525
  • 24
  • 209
  • 285
  • I like the additions. The interesting thing about the `min(which())` method—which I quite like—is that it throws a warning if that simulation does not ever hit a 7, which happens some (I simulated how many times it occurs in 1000 repetitions of the original `replicate()` call, and got an average of 13.9 per 1000.) – JDB Apr 07 '16 at 18:19
  • Yes, `min` will return `Inf` if the condition is never fulfilled. You'll need to increase `n.sim` to ensure that the chain (almost) always reaches 7. I tried with `n.sim=300` and every chain reached 7 in 10,000 simulations. – eipi10 Apr 07 '16 at 18:51
  • Thanks a lot @eipi10 and @JDB! Now I understand the use of loops and `replicate()` a lot better. This is how I solved the problem to avoid `Inf` when computing the average steps: `y<-replicate(1000, min(which(simulation.mc(1,P,100)==7))); y2<-y[is.finite(y)] ; mean(y2)` – Maruska Apr 08 '16 at 10:17
0

It seems like you are trying to create a 100 x 100 data frame, but you're calling an already-created vector instead.

This is why your call of X[i]=rep(0,n.sim) doesn't work. (I'm thinking that you may have meant X[1]=rep(0,n.sim), since you only define i within the loop. You can't fill a column of a vector in the same way that you can a data.frame...

Try:

X <- data.frame(matrix(nrow=100, ncol=100))
for (i in 1:100)
   { X[i]=simulation.mc(1,P,100)
}
JDB
  • 110
  • 6