y = sample(c(rep(1,15), rep(0,35)))
Ntotal=length(y)
datalist = list(y=y,Ntotal=Ntotal)
source("http://www.matsuka.info/univ/course_folder/HDI_revised.txt")

library(rjags)
txt = "
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern( theta )
  }
  theta ~ dbeta( 1, 1 )
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",
                       data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

HDI.plot(mcmcMat)
traceplot(codaSamples)
autocorr.plot(codaSamples,type='l')
gelman.plot(codaSamples)

y1 = sample(c(rep(1,6), rep(0,2)))
y2 = sample(c(rep(1,2), rep(0,5))) 
y = c(y1,y2)
s = c(rep(1,8),rep(2,7))
Ntotal=length(y)
datalist = list(y = y, Ntotal = Ntotal, s = s, Nsubj = 2)

txt = "
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern( theta[s[i_data]] )
}
for ( i_s in 1:Nsubj) {
theta[i_s] ~ dbeta( 2, 2 )
}
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",
                       data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
# checking MCMC
HDI.plot(mcmcMat[,2])
traceplot(codaSamples)
autocorr.plot(codaSamples)
gelman.plot(codaSamples)

x.temp = seq(0,150,0.1)
g1 = dgamma(x.temp, shape =0.01, rate =0.01)
plot(x.temp,g1,type='l')

g2 = dgamma(x.temp, shape =1.56, rate = 0.0312)
plot(x.temp,g2,type='l')

g3 = dgamma(x.temp, shape =6.25, rate = 0.125)
plot(x.temp,g2,type='l')

txt = "
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern( theta[s[i_data]] )
}
for ( i_s in 1:Nsubj) {
theta[i_s] ~ dbeta( omega*(kappa-2)+1,(1-omega)*(kappa-2)+1)
}
omega ~ dbeta(1,1)
kappa <- kappaMinusTwo + 2
kappaMinusTwo ~ dgamma(0.01, 0.01)
}"
writeLines(txt, "~/model.txt")


dat<-read.csv("http://www.matsuka.info/data_folder/TherapeuticTouchData.csv")
y=dat$y
s=as.numeric(dat$s)
Ntotal=length(dat$y)
Nsubj=length(unique(s))
datalist = list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)

jagsModel = jags.model(file="~/model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta","omega","kappa"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

dat<-read.csv("http://www.matsuka.info/data_folder/BattingAverage.csv")
z = dat$Hits
N = dat$AtBats
s = dat$PlayerNumber
c = dat$PriPosNumber
Ncat = 9
Nsubj = 948

datalist = list(z = z, N=N, s=s, c=c, Ncat=Ncat, Nsubj =Nsubj )
txt = "
model {
 for (i_s in 1:Nsubj) {
   z[i_s] ~ dbin( theta[i_s], N[i_s])
   theta[i_s] ~ dbeta(omega[c[i_s]]*(kappa[c[i_s]]-2)+1,
      (1-omega[c[i_s]])*(kappa[c[i_s]]-2)+1)
 }
 for (i_c in 1:Ncat) {
   omega[i_c] ~ dbeta(omegaO*(kappaO-2)+1,(1-omegaO)*(kappaO-2)+1)
   kappa[i_c] <- kappaMinusTwo[i_c]+2
   kappaMinusTwo[i_c] ~ dgamma(0.01,0.01)
 }
 omegaO ~ dbeta(1,1)
 kappaO <- kappaMinusTwoO+2
 kappaMinusTwoO ~ dgamma(0.01,0.01)
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta","omega","kappa"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[," omega[1]"]) #pitcher
HDI.plot(mcmcMat[," omega[2]"]) #catcher
HDI.plot(mcmcMat[," omega[1]"]) #1st base