# JAGS and “simple” GIBBS sampling
# "simple" GIBBS sampling example # initialization n.iter = 10000; sigma = 0.1; counter = 0 z = c(6, 2); N = c(8, 7) a = c(2, 2); b = c(2, 2); n.par = 2 th.hist = matrix(0, nrow = n.iter*n.par, ncol = n.par) theta = runif(2) # function to calc. prob. move prob.gibbs <- function(theta, proposed, N, z, a, b, idx){ p.th=dbeta(theta[idx], z[idx]+a[idx], N[idx]-z[idx]+b[idx]) p.pro=dbeta(proposed, z[idx]+a[idx], N[idx]-z[idx]+b[idx]) return(p.pro/p.th) } # main loop for (i.iter in 1:n.iter){ for (i.par in 1:n.par){ proposed = theta[i.par] + rnorm(1, mean=0, sd=sigma) if (proposed > 1) {proposed = 1} if (proposed < 0) {proposed = 0} p.move= min(1, prob.gibbs(theta, proposed, N, z, a, b, i.par)) if (runif(1) < p.move){ theta[i.par] = proposed } counter = counter + 1 th.hist[counter, ] = theta } } par(mfrow=c(3,1)) HDI.plot(th.hist[,1]) HDI.plot(th.hist[,2]) plot(th.hist[,1],th.hist[,2], type='o',cex=0.1,xlim = c(0,1),ylim=c(0,1)) par(mfrow=c(1,1)) # with JAGS # model model { for (i in 1:Ntotal) { y[i] ~ dbern(theta[ coin[ i ] ]) } for (c in 1:Ncoin) { theta[ c ] ~ dbeta(2, 2) } } # R command y = c(rep(1,6), rep(0,2), rep(1,2), rep(0,5)) coin = c(rep(1,8), rep(2,7)) dataList = list(Ntotal=length(y), y = y, coin = coin, Ncoin = 2) parameters = c("theta") model = jags.model( "~/research/oka/DBDA/ch07/twoVarMCMC.tex", data = dataList, n.chains = 3, n.adapt = 1000) update(model, n.iter=500) CS = coda.samples(model, variable.names = parameters, n.iter = 10000, thin = 5) mcmcMat<-as.matrix(CS) par(mfrow = c(3, 1)) plot(mcmcMat[,1], mcmcMat[,2], type='o',ylim = c(0,1), xlim = c(0,1)) HDI.plot(mcmcMat[,1]) HDI.plot(mcmcMat[,2])