样本 [i-1] 没有使用 Metropolis 填充马尔可夫链蒙特卡洛的 for 循环？答案

【问题标题】：samples[i-1] not populating in for loop for Markov-Chain Monte Carlo using Metropolis?样本 [i-1] 没有使用 Metropolis 填充马尔可夫链蒙特卡洛的 for 循环？
【发布时间】：2020-01-15 18:51:18
【问题描述】：

我正在尝试基于 Metropolis 算法在 R 中为马尔可夫链蒙特卡罗采样生成一个函数。该函数需要接受目标密度函数 (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)

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

【问题讨论】：

samples 的长度为 1，应该是 samples.n + 1。尝试将MCMC.sample 正文的第一行替换为samples <- numeric(length=samples.n + 1)。
非常感谢！！！这是救命稻草。我最终只使用了samples <- numeric(length=samples.n) 而不是length=samples.n +1，因为如果我使用参数samples.n = 1000，添加+1 会得到1001 的输出。
没错，但第一个不会是样本，对吗？或者，您可能会将其视为老化。

标签： r for-loop montecarlo markov-chains

【解决方案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")

【讨论】：

【解决方案2】：

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

【讨论】：