#  BRugs code to compute standardized residuals and CPO's for "stacks" data
#   written by Haijun Ma, Fall 2006
#   modifed by BPC, 2/9/07

# NOTE: Don't forget to change line 2 below to your own working directory,
#   *OR* use "File-Change dir" in R to specify the directory containing
#   your copies of the BRugs and WinBUGS files!

# ALSO:  *DON'T* save these files by right-clicking in Explorer; instead, open 
#   the file and do "File-Save As", *OR* (better yet) open the file and copy-paste 
#   it into a .txt file to ensure proper decoding of carriage returns !!!


oldwd<-getwd() #save your old working directory
setwd("K:\\book\\CL3\\BRugs") #change your working directory


Y<-c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15)
N<-length(Y); p<-3

for(i in 1:N){Ytemp<-Y; Ytemp[i]<-NA; dput(pairlist(p=p,N=N,J=i,
Y=Ytemp, Ytrue=Y), paste("data",i,".txt",sep=""), control=NULL)}
# NOTE: You can also use the "bugsData" function here

mysresid_num<-rep(0,N); myEy2<-rep(0,N); myCPO<-rep(0,N); mymu<-rep(0,N)

library(BRugs)

system.time(     # this is just to time the outer looping procedure
for(i in 1:N){

modelCheck("model.txt")

modelData( paste("data",i,".txt",sep=""))
modelData( paste("dataX.txt",sep=""))
modelCompile(numChains=2)
modelInits(c("Inits1.txt","Inits2.txt"))
modelGenInits()
# NOTE:  You can also use the "bugsInits" function here

modelUpdate(5000)
samplesSet(c("beta0","beta","sigma","tau")) 
summarySet(c("sresid_num","Ey2","CPO","mui")) #running mean and quantiles are calculated. Samples are not stored.
modelUpdate(5000)

#samplesHistory("*", mfrow = c(4, 2))
#samplesBgr("*")
#samplesStats(c("beta0","beta","sigma","tau"))
#dput(samplesStats(c("beta0","beta","sigma","tau")), paste("stats",i,".txt",sep=""))

#record the posterior means of the stats. 
mysresid_num[i]<-summaryStats("sresid_num")[1,1]
myEy2[i]<-summaryStats("Ey2")[1,1]
myCPO[i]<-summaryStats("CPO")[1,1]
mymu[i]<-summaryStats("mui")[1,1]

samplesClear(c("beta0","beta","sigma","tau"))
summaryClear(c("sresid_num","Ey2","CPO","mui"))
}   #  end of outer loop
)   #  end of timing function

mysresid_num
myEy2
myCPO
mymu

setwd(oldwd) #set working directory back to the original one

### Results:
mysresid<-mysresid_num/sqrt(myEy2 - mymu^2)
mysresid

#> mysresid
# [1]  1.08590397 -0.62871474  1.47531084  1.82059477 -0.48356492 -0.87115320
# [7] -0.75114767 -0.41827115 -0.96594276  0.39152783  0.79760128  0.86433071
#[13] -0.42486725 -0.01430217  0.70519663  0.26155406 -0.54451065 -0.13737798
#[19] -0.17813944  0.40264500 -3.02584717

#> myCPO
# [1] 0.122700 0.188400 0.082620 0.047180 0.244200 0.179400 0.184900 0.229700
# [9] 0.158700 0.233700 0.183300 0.164600 0.236400 0.254300 0.192000 0.255500
#[17] 0.185500 0.259300 0.255000 0.249800 0.004631

#  might even draw plots of these:

par(mfrow=c(2,1))
plot(mysresid);  abline(h=0)
plot(myCPO)
mtext("Exact standardized residuals and CPOs, stack loss data",lty=2,line=-2,side=3,outer=T)

#  end of BRugs program