如何用MHadaptive的MCMC拟合多元正态分布解决协方差非正定问题

• 现有代码直接估计协方差矩阵全部9个元素，存在冗余参数，且MH采样的提议分布没有正定约束，一旦生成的矩阵不符合正定要求，dmvnorm和逆威沙特密度函数会直接返回错误导致程序崩溃
• 此前尝试的Cholesky参数化失败是参数约束、先验设置错误导致，方法本身是协方差矩阵估计的标准最优方案

采用Cholesky上三角参数化，仅估计6个独立的Cholesky因子参数+3个均值参数，从根源保证协方差矩阵恒正定：

library(mvtnorm)
library(MHadaptive)
library(LaplacesDemon)

# 生成模拟数据
x <- rmvnorm(100, c(20,65,-44), diag(3))

# 对数后验函数
logpost <- function(parvect, data){
  # 参数拆分：前3个是Cholesky上三角对角线的对数值，后3个是上三角非对角线元素，最后3个是均值
  chol_diag <- exp(parvect[1:3]) # 指数变换保证对角线恒正
  chol_offdiag <- parvect[4:6]
  mu <- parvect[7:9]
  
  # 重构上三角Cholesky因子
  U <- matrix(0, 3, 3)
  diag(U) <- chol_diag
  U[upper.tri(U)] <- chol_offdiag
  sigma <- crossprod(U) # 等价于t(U)%*%U，必然正定
  
  # 先验：协方差用逆威沙特先验，均值用宽正态先验
  log_prior_sigma <- dinvwishart(sigma, nu = 3, S = diag(3), log = TRUE)
  log_prior_mu <- sum(dnorm(mu, 0, 1e5, log = TRUE))
  log_prior <- log_prior_sigma + log_prior_mu
  if(is.infinite(log_prior)) return(-Inf)
  
  # 似然计算
  log_lik <- sum(dmvnorm(data, mean = mu, sigma = sigma, log = TRUE))
  if(is.infinite(log_lik)) return(-Inf)
  
  return(log_lik + log_prior)
}

# 初始值设置：用样本协方差的Cholesky分解作为初始值，保证初始点合法
init_chol <- chol(cov(x))
init_par <- c(log(diag(init_chol)), init_chol[upper.tri(init_chol)], colMeans(x))
par_names <- c(paste0("chol_diag", 1:3), paste0("chol_off", 1:3), paste0("mu", 1:3))

# 运行MH采样
fit <- Metro_Hastings(li_func = logpost, pars = init_par,
                      par_names = par_names, data = x)

# 后处理：采样完成后把Cholesky参数转换回协方差矩阵，即可得到最终参数估计

• 删除了协方差矩阵的冗余参数：3维对称正定矩阵仅需要6个独立参数，不需要估计全部9个元素，避免参数冗余导致的采样效率低下
• Cholesky对角线元素做指数变换，保证重构的协方差矩阵永远正定，完全杜绝非正定报错
• 初始值直接使用样本统计量的Cholesky分解结果，避免初始点非法导致的采样启动失败

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