#########################
#generate a data set in R
#########################
library(truncnorm)
library(mnormt)

dr<-paste('F:/test/',sep='')
setwd(dr)

I=150
J=100

a=rnorm(J,0,sd=.5)
b=rtruncnorm(J,a=0,mean=1.2,sd=.3)

u=rep(NA,I)
for(i in 1:I){
t1=runif(1)
if (t1<.5) u[i]=rnorm(1,-1,.5)
else u[i]=rnorm(1,1,.5)}
u=u/sqrt(1.25) 

d=array(rep(0,I),dim=c(I,1))
for(i in 1:I){
if (runif(1)<pnorm(u[i])) {d[i]=1}
}

w=array(rep(0,I*J),dim=c(I,J))
for(i in 1:I){
  for (j in 1:J){
  if (runif(1)<pnorm(a[j]+b[j]*u[i])) {w[i,j]=1}
  }
}

wdu=cbind(w,d) 
write.table(wdu,"wdu.txt",row.names=FALSE,col.names=FALSE)


######################
#R codes to call JAGS
#######################
library(rjags)

wdu=read.table("wdu.txt")
m=nrow(wdu)
n=ncol(wdu)-1
w=wdu[,1:n]
d=wdu[,n+1]


forJags<-list(D=d,W=w,m=m,n=n)
niter=5000
#initial values
ini.u=rnorm(m)
ini.a=rnorm(n)
ini.b=rexp(n,rate=1)
inits<-list(u=ini.u,a=ini.a,b=ini.b,mua=0,mub=0,taua=.1,taub=.1)

#call JAGS 
foo<-jags.model(file="model.bug",data=forJags,inits=inits,n.chains=2)
out<-coda.samples(model=foo,variable.names=c("u","a","b","mua","taua","mub","taub","wpred"),n.iter=niter,thin=1)
dic.pd=dic.samples(model=foo,n.iter=niter,thin=1,"pD")
dic.pd

#summarize output
iburn=1000
respred=out[[1]]
para=respred[,1:(m+2*n+4)]
meansum=colMeans(para[iburn:niter,])
round(c(meansum[c(1:4,(n+1):(n+4),(2*n+5):(2*n+8),(2*n+1):(2*n+2))],1/sqrt(meansum[(2*n+3):(2*n+4)])),3)

ind=4001
pred=respred[ind,(m+2*n+5):ncol(respred)] 
Wp=matrix(pred,nrow=m,ncol=n)
table(Wp[d==0,])
table(Wp[d==1,])
w=as.matrix(w)
table(w[d==0,])
table(w[d==1,])

###################################################################
#the below codes are saved in a separate file named as "model.bug"
##################################################################
model{
 for (i in 1:m){
  D[i]~ dbern(p[i])
  probit(p[i]) <- u[i]
  u[i]~ dnorm(0,1)

  for (j in 1:n){
  probit(q[i,j]) <- a[j]+b[j]*u[i]  
  W[i,j]~ dbern(q[i,j])
  wpred[i,j]~dbern(q[i,j])
  } 
 }
 
 for (j in 1:n) {
 b[j]~dnorm(mub,taub)T(0,)
 a[j]~dnorm(mua,taua)}

 mua~dnorm(0,.001)
 mub~dnorm(0,.001)
 taua~dgamma(.001,.001)
 taub~dgamma(.001,.001)

}