#############################################
# Example p84
#############################################

   # likelihood in Ex p82
logistic <- function(th)
{
  exp(th)/(1+exp(th))
}

likl.82 <- function(a,b)
{
  x <- c(-0.863,-0.296,-0.053,0.727)
  n <- c(5,5,5,5)
  y <- c(0,1,3,5)
  
  f <- 1
  for(i in 1:4){
    lxi <- logistic(a+b*x[i])
    f <- f * lxi^y[i] * (1-lxi)^(n[i]-y[i])
  }
  f
}

   # posterior in Example p82
post.82 <- function(a,b)
{
  likl.82(a,b)
}

  
  # simulating a sample from the posterior
rpost.82 <- function(n)
{
  i <- sample(1:length(pc),size=n,replace=T,prob=pc)
  ia <- c(row(p))[i]
  ib <- c(col(p))[i]
  a <- agrid[ia]
  b <- bgrid[ib]
  return(cbind(a,b))        # return a n x 2 matrix of a,b vectors
}
  
ld95 <- function(th)
{
  if (is.matrix(th)){ # if passed a whole matrix with 1st col = alpha, 
		      # 2nd col = beta
    alpha <- th[,1]
    beta <- th[,2]
  }
  else{               # if passed only one pair (alpha,beta)
    alpha <- th[1]
    beta <- th[2]
  }
  lambda <- (log(0.95/0.05)-alpha)/beta
  lambda
}

#############################################
# Main commands
#############################################

agrid <- seq(-5,10,length=50)
bgrid <- seq(-10,40,length=50)
p <- matrix(0,nrow=50,ncol=50)    # initialize matrix for posterior
                                    # eval's over agrid x bgrid

# note that post.82 can take a whole vector at a time..
# use that only if you clearly understand what's happ'g
for (i in 1:50){
  p[i,] <- post.82(agrid[i],bgrid)
}
pc <- c(p)  # all cell p's in one vector
image(agrid,bgrid,p,xlab="ALPHA",ylab="BETA") 
contour(agrid,bgrid,p,labex=0,xlab="ALPHA",ylab="BETA")

# marginal in alpha
p.agrid <- apply(p,1,sum)  

# generate a posterior sample
th <- rpost.82(500) # simulate 200 draws from posterior
plot(th[,1],th[,2],xlab="ALPHA",ylab="BETA",xlim=c(-5,10),ylim=c(-10,40),
     pch="o")
contour(agrid,bgrid,p,add=T,labex=0)

# posterior of LD95
lambda <- ld95(th)
hist(lambda,col=0,xlab="LD95",ylab="P(LD95|x)",prob=T)