Metropolis Hastings MCMC with R

I am trying to implement a simple MCMC using the MH algorithm with the R problem, is that I get this error (I was trying to compute the alpha and it is not an NA problem)

Error in if (runif(1) <= alpha) { : missing value where TRUE/FALSE needed

      

here is my function can anyone spot the problem?

    PoissonMetropolisHastingRW = function(n=100000,lambda=10,p=0.5,x0=0){

  x=rep(0,n); y=0; alpha = 0

  x[1]=x0
  for(i in 2:n){

    if (x[i-1] == 0){
      y = sample(c(0,1),1, prob=c(0.5,0.5))
      alpha = min(1,((lambda^y)*x[i-1]*p)/((lambda^x[i-1])*y*(1-p)))
      #alpha = min(1, ( ((lambda^y)*x[i-1])/( (lambda^x[i-1])*y) )*(p/(1-p)) ))
      if(runif(1)<=alpha) {x[i]=y}
      else {x[i]= x[i-1]}

    }
    if (x[i-1] > 0){
      y = sample(c(x[i-1]-1,x[i-1]+1), 1, prob=c(1-p,p))
      alpha = min(1,((lambda^y)*x[i-1]*p)/((lambda^x[i-1])*y*(1-p)))
      #alpha = min(1, (((lambda^y)*x[i-1]/((lambda^x[i-1])*y))*(p/(1-p))))
      if(runif(1) <= alpha) {x[i]=y}
      else {x[i]= x[i-1]}
    }
  } 
  return(x)



}

      

+3


source to share


1 answer


If it y

turns out to be 0 (and with a 0.5 probability for each iteration it will happen with confidence), then it alpha

is 0/0 (because x[i-1] == 0

). He gives you NaN

. The condition something <= NaN

provides NA

.



+3


source







All Articles