Search code examples
rfor-loopnested-loopsmarkov-chainsmarkov

For Loop error in r


I am trying to simulate a markov chain by using code that wasnt made for it. My for loop is not correct and I continue to get error messages but I kept track of the brackets and syntax but I cannot quite put my finger on the issue. The matrix for the markov chain is markov. I am trying to plot each point by appending to the vector a but I dont think my loops are doing exactly what I want.

I am using the random number generator in order to test which state to move to and from. I understand there are easier ways to go about writing this code but I would like to try to use this basic foundation in order to get the code right. If you could please help me straighten the code up and figure out which way to nest the loop properly I would be grateful.

x is my starting point The markov chain probabilities are uniformly distributed

x1 <- runif(100, 1.0, 100)
markov <- matrix(c(.5,.3,.2,.4,.5,.1,.2,.4,.4),nrow=3)
m1 <- markov[1:3] 
m2 <- markov[4:6]
m3 <- markov[7:9]
y <- 2
x <- 1
a <- c(1)

while (y <= 100)
{
if(x==1) {
if(x1[y] < m1[1] * 100) {
x <- 1
} else if (x1[y] <= m1[2]*100 + m1[1] * 100 && x1[y] > m1[1] * 100) {
x <- 2
} else {
x <- 3
} 
a <- union(x,c(y))
} else if (x==2) {
if(x1[y] < m2[1] * 100) {
x <- 1
} else if (x1[y] <= m2[2]*100 + m2[1] * 100 && x1[y] > m2[1] * 100) {
x <- 2
} else {
x <- 3
} 
a <- union(x,c(y))
} else {
if(x1[y] < m3[1] * 100) {
x <- 1
} else if (x1[y] <= m3[2]*100 + m3[1] * 100 && x1[y] > m3[1] * 100) {
x <- 2
} else {
x <- 3
} 
a <- union(x,c(y))
y <- y + 1
}
}
plot(c(1:100),a)

Solution

  • Your y <- y+1 is in one too many curly brackets. It is not being called most of the time.

    I think this is what you want:

    x1 <- runif(100, 1.0, 100)
    markov <- matrix(c(.5,.3,.2,.4,.5,.1,.2,.4,.4),nrow=3)
    m1 <- markov[1:3] 
    m2 <- markov[4:6]
    m3 <- markov[7:9]
    y <- 2
    x <- 1
    a <- c(1)
    aa <- c()
    
    while (y <= 100)
    {
      if(x==1) {
        if(x1[y] < m1[1] * 100) {
          x <- 1
        } else if (x1[y] <= m1[2]*100 + m1[1] * 100 && x1[y] > m1[1] * 100) {
          x <- 2
        } else {
          x <- 3
        } 
        a <- union(x,c(y))
      } else if (x==2) {
        if(x1[y] < m2[1] * 100) {
          x <- 1
        } else if (x1[y] <= m2[2]*100 + m2[1] * 100 && x1[y] > m2[1] * 100) {
          x <- 2
        } else {
          x <- 3
        } 
        a <- union(x,c(y))
      } else {
        if(x1[y] < m3[1] * 100) {
          x <- 1
        } else if (x1[y] <= m3[2]*100 + m3[1] * 100 && x1[y] > m3[1] * 100) {
          x <- 2
        } else {
          x <- 3
        } 
        a <- union(x,c(y))
      }
      aa <- c(aa,x)
    
      y <- y + 1
    }
    plot(aa)