R语言中MCMC估计正态分布参数不收敛的优化建议求助

现有代码的核心问题 #

• 致命语法错误：对数似然计算逻辑完全错误，你写的sum(dnorm(x,mu.old,sd = sigma.old),log = T)是把普通密度求和后给sum传了无用的log参数，正确写法是把log=T放到dnorm的参数里，先算每个点的对数似然再求和。
• 缺失先验项：Metropolis接受比只算了似然差，完全没加入你设定的μ和σ²先验的对数密度，等于没有用到先验信息，不符合贝叶斯推断的逻辑。
• 提议步长不合理：μ和σ的提议标准差都设为10，而你的观测数据均值仅0.7左右，σ估计值也不到1，步长过大导致几乎跳不到合理参数区域，接受率会低到可以忽略。
• 多余约束错误：你额外加了mu.new <=0就拒绝的规则，但你设定的μ先验是均值为0的正态分布，允许取负值，这个约束完全不符合你的先验设定。
• 提议方式效率低：直接对σ做随机游走，靠判断正负拒绝的方式效率极低，σ本身有正约束，更适合对log(σ)做随机游走，自动满足正约束不需要额外判断。
• 存储效率低：循环里用rbind拼接结果矩阵，10万次迭代会大幅拖慢运行速度。

现有代码优化方案 #

优化后的代码如下：

library(MCMCpack)
x = c(0.7669,1.6709,0.8944, 1.0321, 0.0793, 0.1033, 1.2709, 0.7798, 0.6483, 0.3256)
# 先验参数
prior_mu_mean = 0
prior_mu_sd = 100
alpha = 5/2
beta = 10/2

# 初始化参数
mu.old = rnorm(1, prior_mu_mean, prior_mu_sd)
sigma.old = sqrt(rinvgamma(1, shape = alpha, scale = beta))
# 计算初始对数后验 = 对数似然 + 对数先验
loglik.old = sum(dnorm(x, mu.old, sd = sigma.old, log = T))
logprior_mu.old = dnorm(mu.old, prior_mu_mean, prior_mu_sd, log = T)
logprior_sigma2.old = dinvgamma(sigma.old^2, shape = alpha, scale = beta, log = T)
logpost.old = loglik.old + logprior_mu.old + logprior_sigma2.old

# 迭代配置
burn.in = 50000
n_iter = 100000
acc = 0
# 提前初始化结果矩阵
res = matrix(NA, nrow = n_iter, ncol = 3)
colnames(res) = c("mu", "sigma", "logpost")

for (i in 1:n_iter){
  # 调整提议步长，对log(sigma)做游走自动满足正约束
  mu.new = mu.old + rnorm(1, 0, 0.2)
  log_sigma.new = log(sigma.old) + rnorm(1, 0, 0.1)
  sigma.new = exp(log_sigma.new)
  
  # 计算新参数的对数后验
  loglik.new = sum(dnorm(x, mu.new, sd = sigma.new, log = T))
  logprior_mu.new = dnorm(mu.new, prior_mu_mean, prior_mu_sd, log = T)
  logprior_sigma2.new = dinvgamma(sigma.new^2, shape = alpha, scale = beta, log = T)
  logpost.new = loglik.new + logprior_mu.new + logprior_sigma2.new
  
  # Metropolis接受判断，减去log(sigma)的雅可比修正项
  r = exp(logpost.new - logpost.old + log(sigma.new) - log(sigma.old))
  u = runif(1)
  if (u < r){
    mu.old = mu.new
    sigma.old = sigma.new
    logpost.old = logpost.new
    acc = acc + 1
  }
  res[i,] = c(mu.old, sigma.old, logpost.old)
  
  # 每1万次迭代打印接受率
  if (i %% 10000 == 0) {
    cat("迭代", i, "当前接受率：", round(acc/i*100,2), "%\n")
    plot(res[1:i,1], type = 'l', main = "mu trace plot")
  }
}
# 去掉burn-in再做后续推断
res_post = res[-(1:burn.in),]

更高效的MCMC实现方法 #

• 吉布斯采样：这个正态模型属于共轭贝叶斯模型，可以直接推导出μ和σ²的完全条件分布，不需要调步长，接受率100%，收敛速度远快于随机游走Metropolis。μ的完全条件是正态分布，σ²的完全条件是逆伽马分布，直接迭代采样即可，不会出现收敛问题。
• 调用成熟包实现：可以直接用MCMCpack包的MCMCregress函数，或者rstan包拟合模型，不需要自己手写采样逻辑，自动处理先验设置和收敛诊断，稳定性更高。

内容的提问来源于stack exchange，提问作者Song Calderone Zhang