## df.R computes the posterior degrees of freedom rho (dataset 1) in Table 2 
## of Lu, Hodges and Carlin paper using both approximate and exact method
## The details on how to call WinBUGS from R are on
## http://www.stat.columbia.edu/~gelman/bugsR/

y1<-c(35,11,17,14,16,16,24,51,45,9,22,21,9,5,26,8,23,15,14,52)
true.delta1<-c(0.6679837,-0.7090093,-0.08922068,-0.2122527,-0.01051623,-0.54689,-0.02289912,0.8273036,0.7088072,-0.4494048,-0.09783185,-0.01503536,-0.6717945,-1.257309,0.346506,-0.8665467,0.6533667,-0.6930958,-0.6315436,1.205493)

data<-list(y1=y1,N=20)
inits1<-list(tau1 = 2, mu1 = 0, delta1=true.delta1)
inits<-list(inits1, inits1, inits1)
parameters <- c("tau1", "lambda1")
sims<- bugs (data, inits, parameters, "simple-poisson-model.txt",DIC=F,digits.summary=4,plot.summary=F,print.summary=F,n.chains=3, n.iter=1000)

coda<-read.table("coda1.txt")[,2]
coda.index<-read.table("codaIndex.txt")

sigma1<-1/coda[(coda.index[21,2]:coda.index[21,3])]
lambda1<-matrix(0, ncol=20, nrow=coda.index[1,3])

for( i in 1:coda.index[1,3]){
    ind<-coda.index[2:21,2]+i-1
    lambda1[i,]<-coda[ind]
}

appx.df<-function(y,sigma){
    phi<-sigma*y
    w<-1/(1+phi)
    rho<-sum(1-w+(1-w)*w/(N*(1-mean(w))))
    return(rho)
    }

exact.df<-function(lambda, sigma){
    sigma.i<-1/lambda
    phi<-sigma/sigma.i
    wi<-1/(1+phi)
    rho<-sum(1-wi+(1-wi)*wi/(20*(1-mean(wi))))
    return(rho)
}       

appx.rho<-exact.rho<-rep(0,length=coda.index[1,3])
for( i in 1:coda.index[1,3]){
    appx.rho[i]<-appx.df(y1,sigma1[i])
    exact.rho[i]<-exact.df(lambda1[i,], sigma1[i])
}

c(mean(appx.rho), quantile(appx.rho, c(0.025,0.5,0.975)))
c(mean(exact.rho),quantile(exact.rho, c(0.025,0.5,0.975)))