モデル定義ファイルの作成

file_name <- "/tmp/TEMPmodel.txt"

model_string <- "
model {
  theta ~ dbeta(1, 1)
  for(i in 1:Ntotal) {
    y[i] ~ dbern(theta)
  }
}
"
writeLines(model_string, con = file_name)

モデルの作成

y <- c(rep(1, 32), rep(0, 10))

jags_model <- rjags::jags.model(
  file = file_name,
  data = list(
    y = y,
    Ntotal = length(y)
  ),
  inits = list(
    theta = sum(y) / length(y)
  ),
  n.chains = 4
)

MCMC サンプルの生成

# Burn-In
update(jags_model, 500)

# 1000x4 サンプルの生成
coda_samples <- rjags::coda.samples(jags_model, variable.names = "theta", n.iter = 1000)

summary(coda_samples)
#Iterations = 501:1500
#Thinning interval = 1 
#Number of chains = 4 
#Sample size per chain = 1000 
#
#1. Empirical mean and standard deviation for each variable,
#   plus standard error of the mean:
#
#          Mean             SD       Naive SE Time-series SE 
#      0.750295       0.064094       0.001013       0.001028 
#
#2. Quantiles for each variable:
#
#  2.5%    25%    50%    75%  97.5% 
#0.6139 0.7095 0.7547 0.7953 0.8641

可視化

plot(coda_samples)