samples[i-1]没有填充到使用Metropolis的马尔可夫链蒙特卡罗的for循环中？

Question

我正在尝试在R中产生一个函数，用于基于Metropolis算法的马尔可夫链蒙特卡罗抽样。该函数需要接受目标密度函数(PDF)、提出下一步的函数、起点和要评估的步数作为参数。输出应该是一个长度等于步数的向量。

问题是我在尝试使用函数In rnorm(1, x, 0.5) : NAs produced时收到多个警告

我认为问题可能出在我试图定义当前步骤的方式：samples[i-1]不返回值。我不知道为什么会这样。我将samples[1]设置为进入函数的起始点，然后对于i in 2:samples.n，samples[i-1]应该返回前一个值，不是吗？

我试过单独使用函数propose，它工作得很好。这让我相信问题出在for循环中函数propose的输入samples[i-1]。

PDF.beta <- function(x) dbeta(x, 12, 6)

propose <- function(x) rnorm(1, x, 0.5)

MCMC.sample <- function(target.PDF, prop.func, startx, samples.n) {
  samples <- numeric()
  samples[1] <- startx
  for(i in 2:samples.n)
    {
    proposed.step <- prop.func(samples[i-1])
    ifelse(runif(1) < target.PDF(proposed.step)/target.PDF(samples[i-1]),
      samples[i] <- proposed.step,
      samples[i] <- samples[i-1])
  }
  return(samples)
}


beta.MCMC <- MCMC.sample(target.PDF = PDF.beta, prop.func = propose, startx = 4, samples.n = 1000)

我希望看到一个长度为samples.n的向量是从函数MCMC.sample中生成的，给定了beta.MCMC中的输入。相反，samples的输出只有一个值: 4，这是我输入的初始值。我还得到了许多错误消息：In rnorm(1, x, 0.5) : NAs produced。

Answer 1

发布一个完整的正确答案。

主要问题是samples的长度不合适，但我最终还是重写了函数的一部分，因为它似乎没有很好地定义当target.PDF(proposed.step)/target.PDF(samples[i-1])返回NaN时应该发生什么

MCMC.sample <- function(target.PDF, prop.func, startx, samples.n) {
  samples <- numeric(length=samples.n)
  samples[1] <- startx
  for(i in 2:(samples.n)) {
    proposed.step <- prop.func(samples[i-1])
    target.r <- runif(1) < target.PDF(proposed.step)/target.PDF(samples[i-1])
    samples[i] <- ifelse(
      is.na(target.r) | target.r,
      proposed.step, 
      samples[i-1])
  }
  samples
}

set.seed(1)
beta.MCMC <- MCMC.sample(PDF.beta, propose, 4, 1000)

plot(beta.MCMC, type="l")

​

​

Answer 2

在我看来，问题是当目标分布商的分母很小时，商很大。处理这种情况的方法是，当例程产生大于1.0的商(即分母非常小)时，使用if语句将商设置为1.0，然后生成随机数。