How to restart iteration in R loop?

LifeisGood94 0
Viewed 65

I made a loop to generate a Markov Chain. If the proposal does not satisfy a condition, I want to restart the iteration with a new proposal? Is there a way to do this? My current code is shown below for reference. Currently, it sets the current chain's value to the previous one. But I don't want that. I want it to just restart the "i". So if i=2, and the condition in line 4 is not satisfied, I then want it to stay at i=2 until it is satisfied. Thanks in advance.

ABC_MCMC<-function(n){
      for (i in 2:n){
        prop<-rnorm(1,mean=chain[i-1],sd=1)
        if (ABC(prop)==T & prop>=0){
          h_ratio<-(dgamma(prop,shape=prior_alpha,rate=prior_beta)/dgamma(chain[i-1],shape=prior_alpha,rate=prior_beta))*
                   (dnorm(x=chain[i-1],mean=prop,sd=1)/dnorm(x=prop,mean=chain[i-1],sd=1))
          u<-runif(1)
          if (min(1,h_ratio)>u) {chain[i]=prop} else {chain[i]=chain[i-1]}
        }
        else{chain[i]=chain[i-1]}
      }
      return(chain<<-chain)
    }
markov-chains r
1 Answers

This is more of a comment than of an answer but to keep the code formatting I'm posting as an answer.

Replace the code inside the for loop for the code below.

while(TRUE) {
  prop <- rnorm(1, mean = chain[i - 1L], sd = 1)
  if (ABC(prop) && prop >= 0) {
    h_ratio<-(dgamma(prop,shape=prior_alpha,rate=prior_beta)/dgamma(chain[i-1],shape=prior_alpha,rate=prior_beta))*
      (dnorm(x=chain[i-1],mean=prop,sd=1)/dnorm(x=prop,mean=chain[i-1],sd=1))
    u<-runif(1)
    if (min(1,h_ratio)>u) {chain[i]=prop} else {chain[i]=chain[i-1]}
    break
  } else {chain[i] <- chain[i-1]}
}

Edit

The function below seems to be what is asked for.

ABC_MCMC <- function(n){
  for (i in 2:n){
    # loops until condition (ABC(prop) & prop >= 0) is met
    while(TRUE) {
      prop <- rnorm(1, mean = chain[i-1], sd = 1)
      if (ABC(prop) & prop >= 0) {
        h_ratio <- (dgamma(prop, shape = prior_alpha, rate = prior_beta)/dgamma(chain[i - 1L], shape = prior_alpha, rate = prior_beta)) *
          (dnorm(chain[i - 1L], prop, 1)/dnorm(prop, chain[i - 1L], 1))
        u <- runif(1)
        if (min(1, h_ratio) > u) {
          chain[i] <- prop
        } else {
          chain[i] <- chain[i - 1L]
        }
        break
      }
    }
  }
  # function return value
  chain
}
Rui Barradas0
Related