####################
#generate a data set
####################
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)

#############################################
#Codes Below are for the Bayesian estimation
############################################
library(msm)
library(HI)
library(mnormt)

#create functions for arms
mub_fun<-function(x,taub,b,mu,tau){
J=length(b)
tt=-taub/2*sum((b-x)^2)-J*log(pnorm(x*sqrt(taub)))-tau/2*(x-mu)^2
return(tt)
} 

mub_ind_fun<-function(x,taub,b,mu,tau)(x>-10)*(x<10)

taub_fun<-function(x,mub,b,alphab,betab){
J=length(b)
tt=-x*(sum((b-mub)^2)/2+betab)-J*log(pnorm(mub*sqrt(x)))+(J/2+alphab-1)*log(x)
return(tt)
}

taub_ind_fun<-function(x,mub,b,alphab,betab)(x>0)*(x<50)

#read data in
wdu=read.table("wdu.txt")
I=nrow(wdu)
J=ncol(wdu)-1 
w=wdu[,1:J]
d=wdu[,J+1]
wd=wdu[,1:(J+1)]


#arrays for saving MCMC chains
nsim=5000
iburn=1000
u=array(dim=c(I,nsim))
a=array(dim=c(J+1,nsim))
b=array(dim=c(J+1,nsim))
a[J+1,]=0
b[J+1,]=1

mua=array(dim=c(1,nsim))
taua=array(dim=c(1,nsim))
mub=array(dim=c(1,nsim))
taub=array(dim=c(1,nsim))

sen.sim=matrix(nrow=J,ncol=nsim)
spe.sim=matrix(nrow=J,ncol=nsim)

z=array(dim=c(I,J+1))

#hyperparameters in priors
tau=10E-2
mu=0
alphaa=10E-2
betaa=10E-2
alphab=10E-2
betab=10E-2


a[1:J,1]=rnorm(J)
b[1:J,1]=rexp(J,rate=1)
u[,1]=rnorm(I)
mua[,1]=0
taua[,1]=rgamma(1,1)
mub[,1]=0
taub[,1]=rgamma(1,1)

#MCMC loops
 #sample z's
for (m in 1:(nsim-1)){
for(i in 1:I){
 for(j in 1:(J+1)){
 if (wd[i,j]==1) {z[i,j]=rtnorm(1,a[j,m]+b[j,m]*u[i,m],1,lower=0)}
 else {z[i,j]=rtnorm(1,a[j,m]+b[j,m]*u[i,m],1,upper=0)}
 }
}
 
 #sample u's
vu.temp=1/sum(b[,m]^2)
u.temp=vu.temp*((z-matrix(rep(a[,m],I),I,J+1,byrow=T))%*%b[,m])
u[,m+1]=rnorm(I,u.temp,sqrt(vu.temp))
u[,m+1]=u[,m+1]-mean(u[,m+1])
u[,m+1]=u[,m+1]/sd(u[,m+1])

 #sample a's
va.temp=1/(I+taua[,m])
a.temp=va.temp*(colSums(z-u[,m+1]%*%t(b[,m]))+taua[,m]*mua[,m])
a[,m+1]=rnorm(J+1,a.temp,sqrt(va.temp))
a[J+1,m+1]=0

 #sample b's
vb.temp=1/(sum(u[,m+1]^2)+taub[,m])
b.temp=vb.temp*(t(u[,m+1])%*%(z-matrix(rep(a[,m],I),I,J+1,byrow=T))+mub[,m]*taub[,m])
b[,m+1]=rtnorm(J+1,b.temp,sqrt(vb.temp),lower=0)
b[J+1,m+1]=1

 #calculate sensetivity and specificity
sen.sim[,m+1]=pnorm(b[1:J,m+1]%*%t(u[,m+1])+a[1:J,m+1])%*%pnorm(u[,m+1])/I/mean(pnorm(u[,m+1]))
spe.sim[,m+1]=(1-pnorm(b[1:J,m+1]%*%t(u[,m+1])+a[1:J,m+1]))%*%(1-pnorm(u[,m+1]))/I/(1-mean(pnorm(u[,m+1])))
 
 #sample mu_a and tau_a
vmua.temp=1/(J*taua[,m]+tau)
mua.temp=vmua.temp*(taua[,m]*sum(a[1:J,m+1])+tau*mu)
mua[,m+1]=rnorm(1,mua.temp,sqrt(vmua.temp))
taua[,m+1]=rgamma(1,alphaa+J/2,betaa+1/2*sum((a[1:J,m+1]-mua[,m+1])^2))

 #sample mu_b and tau_b
mub[,m+1]=arms(mub[,m],mub_fun, mub_ind_fun,1,taub=taub[,m],b=b[1:J,m+1],mu=0,tau=.01)
taub[,m+1]=arms(taub[,m],taub_fun,taub_ind_fun,1,mub=mub[,m+1],b=b[1:J,m+1],alphab=1,betab=1)
}

#Write out MCMC chains
write(t(a),"a.txt",ncol=nsim)
write(t(b),"b.txt",ncol=nsim)
write(t(u),"u.txt",ncol=nsim)

s=rbind(mua,mub,taua,taub)
write(t(s),"muabtauab.txt",ncol=nsim)


#claculate posterior means, standard deviations, and quantiles
amean=rowMeans(a[1:J,iburn:nsim])
bmean=rowMeans(b[1:J,iburn:nsim])
umean=rowMeans(u[1:I,iburn:nsim])
asd=apply(a[1:J,iburn:nsim],1,sd)
bsd=apply(b[1:J,iburn:nsim],1,sd)
usd=apply(u[,iburn:nsim],1,sd)
aCI=apply(a[1:J,iburn:nsim],1,quantile,probs=c(.025,.975))
bCI=apply(b[1:J,iburn:nsim],1,quantile,probs=c(.025,.975))
uCI=apply(u[,iburn:nsim],1,quantile,probs=c(.025,.975))

muam=mean(mua[iburn:nsim])
tauam=mean(taua[iburn:nsim])
sdam=mean(sqrt(1/taua[iburn:nsim]))
mubm=mean(mub[iburn:nsim])
taubm=mean(taub[iburn:nsim])
sdbm=mean(sqrt(1/taub[iburn:nsim]))
temp=c(muam,tauam,sdam,mubm,taubm,sdbm)

sen.est=rowMeans(sen.sim[,iburn:nsim])
spe.est=rowMeans(spe.sim[,iburn:nsim])

sd.sen.sim=apply(sen.sim[,iburn:nsim],1,sd)
sd.spe.sim=apply(spe.sim[,iburn:nsim],1,sd)

CI.sen=apply(sen.sim[,iburn:nsim],1,quantile, probs=c(.025,.975))
CI.spe=apply(spe.sim[,iburn:nsim],1,quantile, probs=c(.025,.975))

#write out posterior means, standard deviations, and quantiles
write(round(amean,4),"amean.txt",ncol=J)
write(round(bmean,4),"bmean.txt",ncol=J)
write(round(umean,4),"umean.txt",ncol=I)
write(round(asd,4),"asd.txt",ncol=J)
write(round(bsd,4),"bsd.txt",ncol=J)
write(round(usd,4),"usd.txt",ncol=I)
write(round(temp,4),"muabtauabmean.txt",ncol=J)
write(aCI[1,],"aCIL.txt",ncol=J)
write(aCI[2,],"aCIU.txt",ncol=J)
write(bCI[1,],"bCIL.txt",ncol=J)
write(bCI[2,],"bCIU.txt",ncol=J)
write(uCI[1,],"uCIL.txt",ncol=I)
write(uCI[2,],"uCIU.txt",ncol=I)

write(round(sen.est,4),"senmean.txt",ncol=J)
write(round(spe.est,4),"spemean.txt",ncol=J)
write(round(sd.sen.sim,4),"sensd.txt",ncol=J)
write(round(sd.spe.sim,4),"spesd.txt",ncol=J)
write(CI.sen[1,],"senCIL.txt",ncol=J)
write(CI.sen[2,],"senCIU.txt",ncol=J)
write(CI.spe[1,],"speCIL.txt",ncol=J)
write(CI.spe[2,],"speCIU.txt",ncol=J)