R语言是否有成熟预构建包可实现自定义函数的MCMC条件分布采样

我目前使用R编程语言处理如下需求：给定多元联合概率分布函数，需要通过MCMC方法对其条件概率分布进行随机采样。
目前我没有找到现成的实现方案，只能手动编写R代码完成需求，代码虽然可以运行但无法确认正确性，且编写流程繁琐、效率较低，因此希望找到成熟的R第三方包来实现该功能。以下是我当前的实现步骤：

Part 1 - 背景假设 #

假设存在一个4维多元正态分布：

该多元正态分布P(X, Y, Z, W)的参数如下：

• 均值向量（4×1）：
  5.0060022 3.4280049 1.4620007 0.2459998
• 方差-协方差矩阵（4×4）：

[1,] 0.15065114 0.13080115 0.02084463 0.01309107
[2,] 0.13080115 0.17604529 0.01603245 0.01221458 
[3,] 0.02084463 0.01603245 0.02808260 0.00601568 
[4,] 0.01309107 0.01221458 0.00601568 0.01042365

Part 2 - 自定义目标分布函数 #

我基于R语言编写了对应该4维多元正态分布的函数，代码如下：

#define constants needed for the multivariate normal
sigma1 <- c(0.15065114 , 0.13080115 ,  0.02084463 , 0.01309107 , 0.13080115 , 0.17604529 ,  0.01603245 , 0.01221458 , 0.02084463 , 0.01603245  , 0.02808260 , 0.00601568 , 0.01309107 , 0.01221458 ,  0.00601568 , 0.01042365)
sigma <- matrix(sigma1, nrow=4, ncol= 4, byrow = TRUE)
sigma_inv <- solve(sigma)
sigma_det <- det(sigma)
denom = sqrt( (2*pi)^4 * sigma_det) 

#actual multivariate function is defined below ("target")
target <- function(x,y,z,w)
{
  x_one = x - 5.0060022
  x_two = y - 3.4280049
  x_three = z - 1.4620007
  x_four = w - 0.2459998
    
  x_t = c(x_one, x_two, x_three, x_four)
  x_t_one <- matrix(x_t, nrow=4, ncol= 1, byrow = TRUE)
  x_t_two =  matrix(x_t, nrow=1, ncol= 4, byrow = TRUE)
  
  # In this part, as it's (x-mu)^T * SIGMA * (x-mu)
  num = exp(-0.5 * x_t_two  %*%  sigma_inv  %*%  x_t_one)
    
  answer_1 = num/denom
  return(answer_1)
}

Part 3 - 手动实现MCMC采样 #

假设我需要从该多元正态分布的条件分布*P(X, Y | Z = 2 , W = 1.3)*中随机采样。我手动实现了Metropolis-Hastings蒙特卡洛采样器完成该需求，首先固定原多元正态分布中Z和W的取值，代码如下：

#fix the definitions of w and z as per P(X, Y | Z = 2 , W = 1.3) 
target <- function(x,y)
{
  x_one = x - 5.0060022
  x_two = y - 3.4280049
  x_three = 2 - 1.4620007
  x_four = 1.3 - 0.2459998
    
  x_t = c(x_one, x_two, x_three, x_four)
  x_t_one <- matrix(x_t, nrow=4, ncol= 1, byrow = TRUE)
  x_t_two =  matrix(x_t, nrow=1, ncol= 4, byrow = TRUE)
  
  num = exp(-0.5 * x_t_two  %*%  sigma_inv  %*%  x_t_one)
    
  answer_1 = num/denom
  return(answer_1)
}

随后运行采样器完成条件分布的随机采样，代码如下：

library(mvtnorm)
x = rep(0,10000)
y = rep(0,10000)

x[1] = 3     #initialize; I've set arbitrarily set this to 3 and 1
y[1] =1
for(i in 2:10000){
  current_x = x[i-1]
  current_y = y[i-1]
 
  new = rmvnorm(n=1, mean=c(current_x,current_y), sigma=diag(2), method="chol")   # generate bivariate random numbers
  proposed_x = new[1]
  proposed_y = new[2]

  A = target(proposed_x,proposed_y)/target(current_x,current_y) 
  if(runif(1)<A){
    x[i] = proposed_x       # accept move with probabily min(1,A)
    y[i] = proposed_y       
   
    } else {
    x[i] = current_x        # otherwise "reject" move, and stay where we are
    y[i] = current_y
 
    }
}

采样完成后，调用x和y即可获得MCMC结果，最终估计值如下：

mean(mcmc_output$x)
[1] 6.281715

mean(mcmc_output$y)
[1] 4.63817

（可选）Part 4 - 结果可视化 #

我对采样结果的分布进行了可视化，密度图绘制代码如下：

mcmc_output = data.frame(x,y)
par(mfrow=c(1,2))
plot(density(mcmc_output$y, main = "Density of Y"))
plot(density(mcmc_output$x, main = "Density of X"))

等高线图绘制代码如下：

library(ggplot2)

ggplot(mcmc_output, aes(x = x, y = y)) +
  geom_density_2d_filled() + 
  ggtitle("Contour Plots of the MCMC Estimates")

#

核心问题 #

请问是否有更简便的标准方案，可基于知名的预构建R包完成自定义函数的MCMC采样？有没有合适的R包可以推荐？ #

以下是3个经过广泛验证、生态成熟的R包，可以替代手动编写MCMC采样器的流程，避免手动实现的正确性风险，大幅提升开发效率：

1. MCMCpack #

这是专门为MCMC采样设计的主流工具包，内置了Metropolis-Hastings采样的封装函数MCMCmetrop1R，你只需要传入自定义的对数目标密度函数、初始值、迭代次数等参数即可，不需要自己编写循环逻辑和接受/拒绝判断。
针对你的场景，只需要把目标函数改成返回对数密度值（可以避免数值下溢问题，稳定性远高于直接返回密度值，且MH采样本身不需要归一化常数），直接调用即可，示例代码如下：

library(MCMCpack)
# 改写为对数目标密度，固定Z=2、W=1.3
log_target <- function(params) {
  x <- params[1]
  y <- params[2]
  x_one = x - 5.0060022
  x_two = y - 3.4280049
  x_three = 2 - 1.4620007
  x_four = 1.3 - 0.2459998
    
  x_t = c(x_one, x_two, x_three, x_four)
  x_t_one <- matrix(x_t, nrow=4, ncol= 1)
  
  # 直接返回对数密度，不需要除以归一化常数
  return(-0.5 * c(t(x_t_one)  %*%  sigma_inv  %*%  x_t_one))
}
# 调用采样函数
mcmc_res <- MCMCmetrop1R(log_target, theta.init = c(3,1), mcmc = 10000, thin = 1, tune = 1)
# 导出结果
mcmc_output <- as.data.frame(mcmc_res)
colnames(mcmc_output) <- c("x", "y")

2. rstan #

如果你的场景后续要扩展更复杂的自定义分布，或者需要更高的采样效率、自动调参、收敛诊断功能，可以用Stan生态的R接口rstan。你只需要在Stan语法中定义目标分布的对数密度，其余的采样逻辑（采用效率远高于基础MH的NUTS采样器）、收敛检验、结果导出都由框架自动完成。

3. tmvtnorm（仅针对多元正态条件采样场景） #

如果你的实际需求就是多元正态分布的条件采样，根本不需要使用MCMC：多元正态的条件分布有解析解，直接用tmvtnorm包的rmvnorm函数指定条件值即可直接采样，速度比MCMC高几个数量级，结果完全精确，不存在采样误差。

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