在R中运行Metropolis Hastings算法时处理rho参数NaN值的问题

问题背景 #

我在R中针对双变量正态分布的相关参数rho运行Metropolis-Hastings算法，该分布均值为(0,0)，协方差矩阵为[[1, rho], [rho, 1]]，观测值为N组(x,y)，仅关注rho参数。

后验分布（似然×先验）的表达式如下：

prod(1 / ((2 * pi * (1 - rho^2))^(1/2)) * (exp(-(x^2 - 2 * rho * x * y + 
     y^2) / (2 * (1 - rho^2))))*(1 / (1 - rho^2)^(3/2)))

报错代码 #

运行以下R代码时，metropolis函数的if语句出现**"missing value where TRUE/FALSE needed"**错误：

# 获取数据集
x <- data$device1
y <- data$device2

# 定义目标函数（似然×先验）
target <- function(rho){prod(1 / ((2 * pi * (1 - rho^2))^(1/2)) *
     (exp(-(x^2 - 2 * rho * x * y + y^2) / 
     (2 * (1 - rho^2))))*(1 / (1 - rho^2)^(3/2)))}

# 提议分布采用均匀分布U(rho^(i-1)-k, rho^(i-1)+k)
metropolis <- function(rho) {
  rho2 <- runif(1, rho - 0.01, rho + 0.01)
  
  target_rho2 <- target(rho2)
  target_rho <- target(rho)
  
  if (runif(1) > target_rho2/target_rho) {
    rho2 <- rho
  }

  return(rho2)
}

N <- 10000
rho <- mat.or.vec(1, N)
rho[1] <- 0.1
count <- 1

for(n in 2:N){
  repeat{
    rho[n] <- metropolis(rho[n-1])
    count <- count + 1
    if(rho[n]!=rho[n-1]){break}
  }
}

排查与尝试修复 #

已确认错误源于target_rho2/target_rho的比值出现NaN（0/0的情况）。尝试修改if语句，添加如下条件，但错误仍然存在：

is.nan(target_rho2) || is.nan(target_rho) || runif(1) >   
    target_rho2/target_rho
# 或
is.finite(target_rho2) && 
    is.finite(target_rho) && runif(1) > target_rho2/target_rho

内容的提问来源于stack exchange，提问作者hog883