#Approximating the joint posterior distribution p(theta,mu,sigma,tau|data) for the hierarchical normal model described in section 9.8, pag 285, Gelman book
#list of functions
##LIST OF THE FUNCTION---------------------------------
D_data()
I_initial.values(D)
I_one.iteration(I,D)
#gibbs.chain_markov(2000,100,2,I)
#-------------------------------------------------
##DATA AND STARTING POINTS-------------------------- 
data_function(){
  J_4
					#data-set
  yy_rbind(c(62,60,63,59,NA,NA,NA,NA),
           c(63,67,71,64,65,66,NA,NA),  #use scan command for bigger data-sets
           c(68,66,71,67,68,68,NA,NA),
           c(56,62,60,61,63,64,63,59))


					#prior hyperparameters
AA1          _ 3                        #tausqr and sigmasqr have a IG dist
AA2          _ 3                        #with mean close to rough estimates a
BB1          _ 26                        # a big variance
BB2          _ 10 
  output_list(J=J,yy=yy,AA1=AA1,BB1=BB1,AA2=AA2,BB2=BB2)
return(output)
}
 
initial.values_function(D){
  yy_D$yy
  J_D$J
  AA1_D$AA1
  AA2_D$AA2
  BB1_D$AA1
  BB2_D$BB1
  theta_apply(yy,1,mean,na.rm=T)  
  sigmasqr_mean(c(var(yy[1,!is.na(yy[1,])]),
		  var(yy[2,!is.na(yy[2,])]),
		  var(yy[3,!is.na(yy[3,])]),
		  var(yy[4,!is.na(yy[4,])])))
  mu_mean(theta)
  tausqr_var(theta)
  
  output_list(theta=theta,mu=mu,sigmasqr=sigmasqr,tausqr=tausqr)
  return(output)
}

#####ONE ITERATION OF THE GIBBS-SAMPLER------------------------------------------  
one.iteration_function(I,D){
  J_D$J  
  yy_D$yy
  AA1_D$AA1
  AA2_D$AA2
  BB1_D$BB1
  BB2_D$BB2
  theta_I$theta
  mu_I$mu
  sigmasqr_I$sigmasqr
  tausqr_I$tausqr
  n_NULL
				#Draw thetaj's

  for(j in 1:J){
    n[j]_length(yy[j,!is.na(yy[j,])])
    vvv   _ (1/tausqr+n[j]/sigmasqr)^(-1)
    mmm   _ (mu/tausqr+sum(yy[j,!is.na(yy[j,])])/sigmasqr)*vvv
    theta[j]_rnorm(1,mmm,sqrt(vvv)) 
  }
  

					#Draw sigmasqr
  num_NULL
  for(j in 1:J){
   yysqr_sum(yy[j,]^2,na.rm=T)
    num[j]_yysqr+n[j]*theta[j]^2-2*sum(yy[j,],na.rm=T)*theta[j]
  }


num_sum(num)  
  AA1_AA1+sum(n)/2
  BB1_BB1+num/2
  sigmasqr   _ 1/(rgamma(1,AA1)/BB1)
					
#Draw mu

  mu_rnorm(1,mean(theta),sqrt(tausqr/J))

 #Draw tausqr      
  AA2_AA2  + J/2			
  BB2_ BB2 + .5*sum(theta-mu)^2
  tausqr_1/(rgamma(1,AA2)/BB2)
  return(theta,sigmasqr,mu,tausqr)
}

#MARKOV CHAIN
markov _ function(NN,burnin,skip,I) {
# I: list of parameter and data value as in the output of gibbs.city
# skip: number of iteration between two recorded iteration (must be >= 0)
# NN total number of recorded iteration.
# total length of the chain is burnin + skip * NN
  theta _ matrix(NA,length(I$theta),NN)
  sigmasqr_rep(NA,NN)
  mu_rep(NA,NN)
  tausqr_rep(NA,NN)
  for (m in 1:burnin) { I _ one.iteration(I,D) }
  for (m in 1:NN) { 
# if (skip > 0) 
    for (s in 1:skip) { I _ one.iteration(I,D) }
    I _ one.iteration(I,D)
    theta[,m]   _ I$theta
    sigmasqr[m] _ I$sigmasqr
    mu[m]       _ I$mu
    tausqr[m]   _ I$tausqr
  }

return(theta,sigmasqr,mu,tausqr)
}


###MONITORING THE CHAIN 
monitor _ function(x)
{
        par(mfrow = c(2, 1))
        plot(1:length(x), x, xlab = "ITERATION", type = "l")
        hist(x)
        par(mfrow = c(1, 1))
}


plotmeans_function(w=0,I=gibbs.chain,cc=2){
if (w == 1) postscript("/home/biostats/fdominic/course/plotmeans.ps")
par(oma=c(0,0,2,0))
boxplot(I$theta[1,],I$theta[2,],I$theta[3,],I$theta[4,],I$mu,
names=c("theta1","theta2","theta3","theta4","mu"),style.bxp = "old",cex=cc)
par(oma=c(0,0,0,0))
par(mfrow=c(1,1))
if (w == 1) dev.off()

}
plotsigma_function(w=0,I=gibbs.chain,cc=2){
if (w == 1) postscript("/home/biostats/fdominic/course/plotsigma.ps")
par(oma=c(0,0,2,0))
hist(sqrt(I$sigmasqr),xlab="sigma",density=-1,yaxt="n",cex=cc)
abline(v=mean(I$sigma),lty=2)
par(oma=c(0,0,0,0))
par(mfrow=c(1,1))
if (w == 1) dev.off()

}
plottau_function(w=0,I=gibbs.chain,cc=2){
if (w == 1) postscript("/home/biostats/fdominic/course/plottau.ps")
par(oma=c(0,0,2,0))
hist(sqrt(I$tausqr),xlab="tau",nclass=20,density=-1,yaxt="n",xlim=c(0,15),cex=cc)
par(oma=c(0,0,0,0))
par(mfrow=c(1,1))
if (w == 1) dev.off()

}