3

Rで分析を行うためにrjagsを使用しています(このブログ投稿に基づいています:http://www.sumsar.net/blog/2013/08/bayesian-estimation-of-correlation/)そして質問があります. update() を使用してモデルを 500 サンプルで更新することと、元のモデル作成で n.adapt パラメータを 500 ではなく 1000 に設定することとの間に違いはありますか?

つまり、次のことを行う場合、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