# MCMC

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