如何使用 MCMC 拟合多元正态分布（MHadaptive）答案

【问题标题】：how to fit multivariate normal distribution with MCMC (MHadaptive)如何使用 MCMC 拟合多元正态分布（MHadaptive）
【发布时间】：2021-07-02 11:11:03
【问题描述】：

我想通过 MCMC 抽样来估计多元正态分布的参数。棘手的部分是协方差矩阵。我知道传统的先验选择是逆愿望先验（或精度矩阵的愿望先验）。但是，我不知道该给谁用。当我尝试估计多元正态分布的参数时，MHadaptive 崩溃，因为 sigma 不是正定的。我还用不同的采样器尝试了它，并最终遇到了同样的问题。以下是我迄今为止尝试过的其他一些事情：

使用带有cholesky 参数化的逆愿望图并估计U 的上三角。然后sigma=t(U)%*%U。这不会导致正定问题，但它也不起作用。

改用精度矩阵，也没有解决问题。

当然，我也试过只估计 Sigma 的上三角并重构 Sigma，使其对称。但是，它仍然不是肯定的。

这是我的代码：

library(mvtnorm)
library(MHadaptive)

#data
x <- rmvnorm(100, c(20,65,-44), diag(3))

#define the log prior
logprior <- function(parvect){
    logprior <- LaplacesDemon::dinvwishart(matrix(parvect[1:9],3,3,byrow=TRUE),3, diag(3), log=TRUE) +
      dnorm(parvect[10], 0, 100000, log=TRUE) + 
      dnorm(parvect[11], 0, 100000, log=TRUE) + 
      dnorm(parvect[12], 0, 100000, log=TRUE) 
    return(logprior)
}

#define the log likelihood
LL <- function(parvect, data){
  LL <- sum(mvtnorm::dmvnorm(data, parvect[10:12], matrix(parvect[1:9],3,3,byrow=TRUE), log=TRUE))

  return(LL) 
}

LL_reg <- function(parvect, data){
  LL <- LL(parvect = parvect, data=data)
  logPrior <- logprior(parvect = parvect)
  LL_reg <- LL+logPrior
  return(LL_reg)
}

#inits etc.
df_params <- data.frame(name = c( "covmat[1,1]","covmat[1,2]", "covmat[1,3]",
                                  "covmat[2,1]","covmat[2,2]","covmat[2,3]",
                                  "covmat[3,1]","covmat[3,2]","covmat[3,3]",
                                  "mean[1]", "mean[2]", "mean[3]"),
                        min = c(rep(-Inf,12)),max = c(rep(Inf,12)),init=c(as.vector(t(cov(x))), colMeans(x)))


Metro_Hastings(li_func = LL_reg, pars = df_params$init,
               par_names = df_params$name, data=x)

【问题讨论】：

标签： bayesian mcmc

【解决方案1】：

既然你知道可能性和先验，为什么不试试 Rstan 来获得估计？定义 cov_matrix 或直接将方差-协方差矩阵放在多正态函数的方差部分会限制它为 p.d.

【讨论】：

好吧，我的目标不是适应分布，而是要了解如何解决一般问题以及为什么会发生这种情况。我知道我可以很容易地用 stan 或 jags 做到这一点。