#Example pag 121
 
 log.post.tarone_function(gamma,delta){
   tarone _ matrix(scan("tarone.txt"),byrow=T,ncol=2)
   deaths _ tarone[,1]
   rats   _ tarone[,2]
   alpha  _ exp(gammadelta)/(1 + exp(gamma))
   beta   _ exp(delta)/(1+exp(gamma))
   ldens  _ 0
   for (i in 1:length(tarone[,1])){
     ldens  _ ldens +
       (lgamma(alpha+beta)+lgamma(alpha+deaths[i])+lgamma(beta+rats[i]-deaths[i]))-
         (lgamma(alpha)+lgamma(beta)+lgamma(alpha+beta+rats[i]))
   } 
   ldens - 5/2*log(alpha+beta)+log(alpha)+log(beta)
 }
 
 plot.joint.post.tarone_function(w=0,ll=200,sim){
   if (w == 1) postscript("/home/biostats/fdominic/course/jposterior.tarone.ps")
   gamma _ seq(-2.3,-1.3,length=ll)
   delta  _ seq(1,5,length=ll)
   contours _ seq(.05,.95,.1)
   logdens _outer(gamma,delta,log.post.tarone) 
  dens_exp(logdens-max(logdens))
  contour(gamma,delta,dens,levels=contours,xlab="log(alpha/beta)",ylab="log(alpha+beta)",labex=0,cex=1)
points(log(sim$aaa/sim$bbb),log(sim$aaa+sim$bbb)) 
 mtext("Posterior density",3,line=1,cex=1.2)
 par(oma=c(0,0,0,0))
  par(mfrow=c(1,1))
  if (w == 1) dev.off()}




grid.value.tarone_function(ll=10){
  gamma _ seq(-2.3,-1.3,length=ll)
  delta  _ seq(1,5,length=ll)
  PP    _ matrix(NA,ll,ll)
  for(i in 1:ll){ 
    for(j in 1:ll){
      PP[i,j]_log.post.tarone(gamma[i],delta[j])	
    }}
  MM _ max(PP)
  PP_exp(PP-MM)
  ccc _ sum(PP)
  PP  _ PP/ccc
  return(PP,gamma,delta)
}


  
        sampling.tarone_function(NN=1000,II,gamma,delta){
	  ggg   _ aaa _ ddd _ bbb _ NULL
	  gamma.mar_apply(II,1,sum)
	  delta.mar_apply(II,2,sum) 
	  tarone _ matrix(scan("tarone.txt"),byrow=T,ncol=2)
	  deaths _ tarone[,1]
	  rats   _ tarone[,2]
	  ttt_matrix(NA,length(rats),NN)
	  for(l in 1:NN){
	    uuu_runif(1,0,1)
	    ggg[l]_max(gamma[cumsum(gamma.mar) < uuu])
	    junk_(1:length(gamma))[gamma==ggg[l]]
	    conditional_II[junk,]/sum(II[junk,])
	    uuu_runif(1,0,1)
	    ddd[l]_max(delta[cumsum(conditional) < uuu])
	    aaa[l]  _ exp(ggg[l]+ddd[l])/(1+exp(ggg[l]))
	    bbb[l]  _ exp(ddd[l])/(1+exp(ggg[l]))
	
    for(i in 1:length(deaths)){
	      ttt[i,l]_rbeta(1,aaa[l]+deaths[i],bbb[l]+rats[i]-deaths[i])
	    }
	  }
return(aaa,bbb,ttt)
      }





make.boxplot_function(w = 0,sim){
 if(w == 1) postscript("/home2/biostats/fdominic/course/boxtheta.ps")     
freq_labels(1:71)
theta_sim$ttt
TT_length(theta[,1])  #number of experiments
boxplot(split(t(theta),rep(1:TT,rep(length(theta[1,]),TT))),style.bxp="old",names=freq,srt=45,cex=.7)
par(oma=c(0,0,0,0))
if(w == 1) dev.off()
}