#read a data set in
data<-read.table("simdata.txt")
L<-data[,1]
R<-data[,2]
icen<-data[,3] #an indicator to show the censoring status
x1<-data[,4] 
z<-data[,5:6]  #two repeated observations with measurement error
n=nrow(data)   #sample size
p=2            #dimension of covariate 

l=rep(0,n)
r=rep(0,n)
r[icen==0]=log(R[icen==0])#left censored
for(i in ((1:n)[icen==0])){l[i]=min(-1,(r[i]-.1))}
l[icen==1]=log(L[icen==1])#interval censored
r[icen==1]=log(R[icen==1])
l[icen==2]=log(L[icen==2])#right censored
for(i in ((1:n)[icen==2])){r[i]=max(9,l[i]+.1)}
l[icen==3]=log(L[icen==3])#exactly observed
r[icen==3]=log(R[icen==3])


#specify hyperparameters for priors
nv=.1
atau=.1
btau=.1
alpM=.1
betM=.1
alptau=.1
beltau=.1
alpx=.1
betx=.1
beta0=c(0,0)
Sig0=matrix(c(1,0,0,1),2,2)

nsim=5000 #number of MCMC iterations
y=array(dim=c(n,nsim))
M=array(dim=c(1,nsim))
bet=array(dim=c(p,nsim))
x=array(dim=c(n,nsim))
taue=array(dim=c(1,nsim))
taux=array(dim=c(1,nsim))
mux=array(dim=c(1,nsim))
nclu=array(dim=c(1,nsim))
t=seq(-5,5,.01)
denW=array(dim=c(length(t),nsim))

#specify initial values for MCMC loops
N=20
muh=rep(NA,N)
tauh=rep(NA,N)
U=rep(NA,n)
s=rep(NA,n)
vh=rep(NA,N)
ph=rep(NA,N)
ind=array(dim=c(N,n))

M[1,1]=1
bet[,1]=beta0
x[,1]=rowMeans(z)
mux[,1]=0
taux[,1]=rgamma(1,alptau,bettau)

muh=rnorm(N,0,sqrt(1/nv))
tauh=rgamma(N,atau,rate=btau)

P=rep(1/N,N)
ind=rmultinom(n,1,P) 
for(i in 1:n){s[i]=seq(1,N)[ind[,i]==1]}
nh=rowSums(ind)


for (m in 1:(nsim-1)){

xcov=cbind(x1,x[,m])
y[icen==3,m]=l[icen==3]
for (i in ((1:n)[icen!=3])){
mu=muh[s[i]]-xcov[i,]%*%bet[,m]
sig=1/sqrt(tauh[s[i]])
y[i,m]=rtnorm(1,mu,sig,lower=l[i],upper=r[i])}

for (i in 1:N){
sig2=1/(tauh[i]*nh[i]+nv)
mu=sig2*( tauh[i]*sum((y[,m]+xcov%*%bet[,m])[s==i]))
muh[i]=rnorm(1,mu,sqrt(sig2))
tauh[i]=rgamma(1,atau+nh[i]/2,rate=btau+1/2*sum(((y[,m]+xcov%*%bet[,m]-muh[i])[s==i])^2))}

tempn=sum(nh)-cumsum(nh)
for (i in 1:N){
vh[i]=rbeta(1,nh[i]+1,M[1,m]+tempn[i])
vh[i]=min(1-1e-6,vh[i]) }

tempv=cumprod(c(1,1-vh))
for(i in 1:N){ph[i]=tempv[i]*vh[i]}

M[1,m+1]=rgamma(1,alpM+N,rate=betM-sum(log(1-vh)))


for(i in 1:n){U[i]=runif(1,max=ph[s[i]])}
ustar=min(U)
NN=N
while (sum(ph)<1-ustar){
muh=c(muh,0)
tauh=c(tauh,0)
vh=c(vh,0)
ph=c(ph,0)
muh[N+1]=rnorm(1,0,sqrt(1/nv))
tauh[N+1]=rgamma(1,atau,rate=btau)
vh[N+1]=rbeta(1,1,M[1,m+1])
vh[N+1]=min(1-1e-10,vh[N+1])
ph[N+1]=ph[N]*vh[N+1]*(1-vh[N])/vh[N]
N=N+1}

ind=array(dim=c(N,n))
for (i in ((1:n)[icen==3])){
ind[,i]=rmultinom(1,1,(dnorm(y[i,m],-xcov[i,]%*%bet[,m]+muh,1/sqrt(tauh)))*(ph>U[i]))}
for (i in ((1:n)[icen!=3])){
ind[,i]=rmultinom(1,1,(pnorm(sqrt(tauh)*(r[i]+xcov[i,]%*%bet[,m]-muh))-pnorm(sqrt(tauh)*(l[i]+xcov[i,]%*%bet[,m]-muh)))*(ph>U[i]))}

for(i in 1:n){s[i]=seq(1,N)[ind[,i]==1]}
nh=rowSums(ind)
nclu[1,m]=sum(nh!=0)

for (i in 1:length(t)){denW[i,m]=sum(ph*dnorm(tt[i],muh,1/sqrt(tauh)))}

Sig=solve(solve(Sig0)+t(tauh[s]*xcov)%*%xcov)
Beta=Sig%*%(solve(Sig0)%*%beta0+t(xcov)%*%(tauh[s]*(-z[,m]+muh[s])))

SigE=eigen(Sig)
SigEV=matrix(c(SigE$values[1],0,0,SigE$values[2]),2,2)
SigH=SigE$vectors%*%sqrt(SigEV)

bet[,m+1]=Beta+SigH%*%rnorm(p)

taue[m]=rgamma(1,alptau+2*n/2,rate=bettau+1/2*sum((z-x[,m])^2))

sigx=1/(tauh[s]*bet[2,(m+1)]^2+2*taue[m]+taux[m])
muxx=sigx*(rowSums(z)*taue[m]+tauh[s]*(muh[s]-y[,m]-xcov[,-2]*bet[-2,m+1])*bet[2,m+1]+taux[m]*mux[m])
for(i in 1:n){x[i,m+1]=rnorm(1,muxx[i],sqrt(sigx[i]))}

sigtemp=1/(n*taux[m]+nv)
mutemp=sigtemp*taux[m]*sum(x[,m+1])
mux[m+1]=rnorm(1,mutemp,sqrt(sigtemp))
taux[m+1]=rgamma(1,alpx+n/2,rate=betx+1/2*sum((x[,m+1]-mux[m+1])^2))
}


iburn=1001
nthin=1

muxmean=mean(mux[seq(iburn,(nsim-1),nthin)])
muxsd=sd(mux[seq(iburn,(nsim-1),nthin)])
muxq=quantile(mux[seq(iburn,(nsim-1),nthin)],c(.025,.975))
tauxmean=mean(taux[seq(iburn,(nsim-1),nthin)])
tauxsd=sd(taux[seq(iburn,(nsim-1),nthin)])
tauxq=quantile(taux[seq(iburn,(nsim-1),nthin)],c(.025,.975))

tauemean=mean(taue[seq(iburn,(nsim-1),nthin)])
tauesd=sd(taue[seq(iburn,(nsim-1),nthin)])
taueq=quantile(taue[seq(iburn,(nsim-1),nthin)],c(.025,.975))


bet1mean=mean(bet[1,seq(iburn,(nsim-1),nthin)])
bet2mean=mean(bet[2,seq(iburn,(nsim-1),nthin)])
bet1sd=sd(bet[1,seq(iburn,(nsim-1),nthin)])
bet2sd=sd(bet[2,seq(iburn,(nsim-1),nthin)])
bet1q=quantile(bet[1,seq(iburn,(nsim-1),nthin)],c(.025,.975))
bet2q=quantile(bet[2,seq(iburn,(nsim-1),nthin)],c(.025,.975))
c(bet1mean-2*bet1sd,bet1mean+2*bet1sd)
c(bet2mean-2*bet2sd,bet2mean+2*bet2sd)


denWmean=rowMeans(denW[,seq(iburn,(nsim-1),nthin)])

plot(t,denWmean)