# 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])