3

我正在使用 rjags 在 R 中进行分析（基于这篇博文：http: //www.sumsar.net/blog/2013/08/bayesian-estimation-of-correlation/）我有一个问题. 使用 update() 更新具有 500 个样本的模型与仅将原始模型创建中的 n.adapt 参数设置为 1000 而不是 500 之间有什么区别吗？

换句话说：当我执行以下操作时，run1 和 run2 之间是否存在重要区别：

require(mvtnorm)
require(rjags)

model_string <- "
  model {
    for(i in 1:n) {
      x[i,1:2] ~ dmnorm(mu[], prec[ , ])
    }

    prec[1:2,1:2] <- inverse(cov[,])
    cov[1,1] <- sigma[1] * sigma[1]
    cov[1,2] <- sigma[1] * sigma[2] * rho
    cov[2,1] <- sigma[1] * sigma[2] * rho
    cov[2,2] <- sigma[2] * sigma[2]

    sigma[1] ~ dunif(0, 1000) 
    sigma[2] ~ dunif(0, 1000)
    rho ~ dunif(-1, 1)
    mu[1] ~ dnorm(0, 0.001)
    mu[2] ~ dnorm(0, 0.001)

    x_rand ~ dmnorm(mu[], prec[ , ])
  }

"

mu <- c(10, 30)
sigma <- c(20, 40)
rho <- -0.7
cov_mat <- rbind(c(     sigma[1]^2       , sigma[1]*sigma[2]*rho ),
                 c( sigma[1]*sigma[2]*rho,      sigma[2]^2       ))
x <- rmvnorm(30, mu, cov_mat)

data_list = list(x = x, n = nrow(x))
inits_list = list(mu = c(mean(x[, 1]), mean(x[, 2])),
                  rho = cor(x[, 1], x[, 2]),
                  sigma = c(sd(x[, 1]), sd(x[, 1])))


jags_model <- jags.model(textConnection(model_string), data = data_list, inits = inits_list,
                         n.adapt = 500, n.chains = 3, quiet = T)
update(jags_model, 500)
mcmc_samples <- coda.samples(jags_model, c("mu", "rho", "sigma", "x_rand"),
                             n.iter = 5000)

run1<-summary(mcmc_samples)

jags_model <- jags.model(textConnection(model_string), data = data_list, inits = inits_list,
                         n.adapt = 1000, n.chains = 3, quiet = T)
mcmc_samples <- coda.samples(jags_model, c("mu", "rho", "sigma", "x_rand"),
                             n.iter = 5000)
run2<-summary(mcmc_samples)
4

1 回答 1

5

该函数update为您的马尔可夫链执行老化迭代。如果您未能运行这些老化迭代，那么您的链将不会收敛到其平稳分布，并且您生成的任何样本实际上都不会来自正确的后验分布。简而言之，您绝对需要更新您的链，并且您还应该对最终样本执行收敛检查（跟踪图/acfs），以确保您已经将链烧掉足够长的时间。

JAGS 执行的自适应迭代不会在迭代中燃烧。从文档中，“在此自适应阶段生成的样本序列不是马尔可夫链，因此可能不能用于模型的后验推断。” 在自适应阶段运行模型可以提高效率，但您仍然需要更新/老化您的链，以便您实际上是从正确的后验分布中采样。

于 2015-05-05T17:21:51.247 回答