deviance<-function(y,subject,mu,sd){
   -2*sum(dnorm(y,mu[subject],sd,log=T))
}

ONEWAYRE<-function(y,subject,tau2=1,n.samples=10000,burn=1000){
   #  Fit the model: 
   #    y[i]~N(mu[subject[i]],prec=tau1)
   #    mu[j]~N(mu.bar,tau2)
   #    mu.bar~flat

   n<-length(y)
   n.subs<-max(subject)

   mu<-rnorm(n.subs)
   mu.bar<-mean(mu)
   tau1<-1

   keepmu<-matrix(0,n.samples,n.subs)
   keeptau1<-keepmu.bar<-rep(0,n.samples)
   dev<-LI<-rep(0,n.samples)
   
   for(i in 1:n.samples){
     tau1<-rgamma(1,n/2+0.01,rate=sum((y-mu[subject])^2)/2+0.01)
     
     for(j in 1:n.subs){
       COV<-1/(tau1*sum(subject==j)+tau2)
       MN<-COV*(tau1*sum(y[subject==j])+tau2*mu.bar) 
       mu[j]<-rnorm(1,MN,sqrt(COV))
     }
    
     mu.bar<-rnorm(1,mean(mu),1/sqrt(tau2*n.subs))  

     keepmu[i,]<-mu
     keeptau1[i]<-tau1
     keepmu.bar[i]<-mu.bar

     #Compute the deviance
     dev[i]<-deviance(y,subject,mu,1/sqrt(tau1))
     
     #Compute Laud and Ibrahim's statistic:
     ynew<-rnorm(n,mu[subject],1/sqrt(tau1))
     LI[i]<-sum((y-ynew)^2)     
   }

   mu.hat<-apply(keepmu[burn:n.samples,],2,mean)
   sd.hat<-mean(1/sqrt(keeptau1[burn:n.samples]))

   dbar<-mean(dev[burn:n.samples])
   dhat<-deviance(y,subject,mu.hat,sd.hat)
   pD<-dbar-dhat
   DIC<-dbar+pD

list(mu=keepmu,mu.bar=keepmu.bar,tau1=keeptau1,dev=dev,
     dbar=dbar,dhat=dhat,pD=pD,DIC=DIC,LI=mean(LI[burn:n.samples]))}
  
n.subs<-20
n.reps<-2
subject<-rep(1:n.subs,n.reps)

n.sims<-100
LI<-DIC<-DBAR<-pD<-matrix(0,n.sims,3)

for(s in 1:n.sims){
  mu<-rnorm(n.subs,3,1)
  y<-rnorm(n.reps*n.subs,mu[subject],1)
  fit1<-ONEWAYRE(y,subject,tau2=exp(10))
  fit2<-ONEWAYRE(y,subject,tau2=1)
  fit3<-ONEWAYRE(y,subject,tau2=exp(-10))
  LI[s,]<-c(fit1$LI,fit2$LI,fit3$LI)
  DBAR[s,]<-c(fit1$dbar,fit2$dbar,fit3$dbar)
  DIC[s,]<-c(fit1$DIC,fit2$DIC,fit3$DIC)
  pD[s,]<-c(fit1$pD,fit2$pD,fit3$pD)
}