我正在尝试使用带有R的MH算法来实现一个简单的MCMC问题,就是我遇到了这个错误(我试图计算alpha值,这不是NA问题)

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


这是我的功能,谁能发现问题?

    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)



}

最佳答案

如果y恰好为0(并且每次迭代的概率为0.5,则可以肯定地发生),则alpha为0/0(因为x[i-1] == 0)。它为您提供NaN。条件something <= NaN提供一个NA

关于r - 带有R的Metropolis-Hastings MCMC,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/29698831/

10-14 10:11
查看更多