#############################################
# Hierarchical model, Rat tumors p129
#############################################

G <- 50                                # set grid sizes

data <- scan("rats.data",what=list(1,1))
				       # read in data from "rats.data"
y <- data[[1]]                         # y is the first data column
nj <- data[[2]]                        # nj the second
N <- length(y)

# evaluates log prior for phi=(u,v) in par's u=log(a/b),v=log(a+b)
# note: might not need a transformation for other examples 
# (do need thought for 5.11) 
lprior <- function(u,v)
{
  b <- exp(v)/(1+exp(u)) # transform to alpha,beta
  a <- b*exp(u)
  log(a)+log(b)-5/2*v
}

# evaluates the marginal posterior in phi=(u,v) (in logs)
lmarg.post <- function(u,v)
{
  b <- exp(v)/(1+exp(u)) # transform to alpha,beta
  a <- b*exp(u)
  f <- 0
  for (i in 1:N){
    f <- f+lgamma(a+y[i])+lgamma(b+nj[i]-y[i])-lgamma(a+b+nj[i])
  }
  f <- f + N*(lgamma(a+b)-lgamma(a)-lgamma(b))
  lp <- f+log(a)+log(b)-5/2*v   # add the log prior
  lp                            # return the log posterior
}

# generates n draws from the marginal posterior p(a,b)
# needs 
# p:    matrix of posterior evaluated at ugrid x vgrid
# pc:   matrix p as one big vector
# both are created in the main program below..
marg.rpost <- function(n)
{
  i <- sample(1:length(pc),size=n,replace=T,prob=pc)
  iu <- c(row(p))[i]
  iv <- c(col(p))[i]
  u <- ugrid[iu]
  v <- vgrid[iv]
  return(cbind(u,v))        # return a n x 2 matrix of a,b vectors
}


#############################################
# Begin main commands
#############################################

ugrid <- seq(-2.3,-1.4,length=G)     # grids for phi=(a,b)
vgrid <- seq(1,5,length=G)
lp <- matrix(0,nrow=G,ncol=G)    # initialize matrix for posterior

#############################################
# evaluate log posterior on the grid
for (i in 1:G){
  lp[i,] <- lmarg.post(ugrid[i],vgrid)
}
# divide posterior by max(exp(lp)) (i.e. substract max(lp) from lp) to avoid 
# numerical problems
lp <- lp-max(lp)
image(ugrid,vgrid,lp,xlab="LOG(A/B)",ylab="LOG(A+B)")  # log posterior

p  <- exp(lp)                    # now convert from logs to natural scale
                                 # and plot the posterior
pc <- c(p)                       # all cell p's in one vector
                                 # needed for "marg.rpost()"
contour(ugrid,vgrid,p,labex=0,xlab="LOG(A/B)",ylab="LOG(A+B)") 
image(ugrid,vgrid,p,xlab="LOG(A/B)",ylab="LOG(A+B)") 

#############################################
# simulate from the posterior on (u,v)

# generate from the posterior on (u,v)
uv <- marg.rpost(1000)
plot(uv,xlim=range(ugrid),ylim=range(vgrid),pch="*",xlab="LOG(A/B)",
     ylab="LOG(A+B)")

#############################################
# generate th[j]'s

u <- uv[,1]
v <- uv[,2]
b <- exp(v)/(1+exp(u)) # transform to alpha,beta
a <- b*exp(u)

thjbar <- rep(0,N)     # initialize list of E(th[j]|y)
thlo   <- rep(0,N)     #            list of lower 0.025 quantiles
thhi   <- rep(0,N)
for(j in 1:N){
  thj <- rbeta(1000,a+y[j],b+nj[j]-y[j]) # simulate 1000 draws for th[j]
  thjbar[j] <- mean(thj)
  thlo[j]   <- sort(thj)[25]
  thhi[j]   <- sort(thj)[975]
}
thhat <- y/nj
plot(thhat,thjbar,pch=15,xlab="Yj/Nj",ylab="95% POST INTERVAL",ylim=c(0,0.4),
     xlim=c(0,0.4))
for(j in 1:N){ # add lines from (thhat_j,thlo_j) to (thhat_j,thhi_j)
  lines(c(thhat[j],thhat[j]),c(thlo[j],thhi[j]))
}
abline(0,1,col=2)