# likelihood in Example p. 82 

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)  
 }
 
 agrid<- seq(-5,10,length=200)
 bgrid <- seq(-10,40,length=200)
 p <- matrix(0,nrow=200,ncol=200) # initialize matrix for posterior 
                                                                # eval's over 
agrid x bgrid 
                                                
 #evaluate posterior on grid 
 for(i in 1:200) {
    for(j in 1:200) { 
     p[i,j] <- post.82(agrid[i],bgrid[j]) 
  }
     }

#alternatively --note that post.82 can take a whole vector t a time.. 
#                          use that only if you clearly understand what's 
happening
for(i in 1:200) {
    p[i,] <- post.82(agrid[i],bgrid) 
  }
  
contour(agrid,bgrid,p,xlab="ALPHA", ylab="BETA")
image(agrid, bgrid,p,xlab="ALPHA",ylab="BETA")


# First create the functions that we will need to generate the 
distributions 

##############################################################
#     Function init.rpost.82  
 # requires: agrid, bgrid, p, p.agrid
 # sets up  marginal cdf for alpha:                 cdfa
 #                 cond cdf for beta   | alpha:            cdfba 
 # note: computing cdfba-inf for p(beta | alpha) would be too cumbersome, 
 # would need it on a whole grid of alphas ! 
 #############################################################
 init.rpost.82 <- function()
 {
 da <- agrid[2]-agrid[1]              #step size on a grid 
 k <- sum(p.agrid*da)        # probability in that interval
 cat("computing cdfa.inv....\n")
 cdfa <<- rep(0,200)   # initialize marginal cdf for alpha
 cdfa[1] <<- p.agrid[1]*da/k
 
 for ( i in 2:200)          #compute marginal cdf for alpha 
   cdfa[i] <<- cdfa[i-1] + da*p.agrid[i]/k
      
      cdfba <<- matrix(0, nrow=200,ncol=200)   # initialize cdf for beta 
given alpha 
         db <- bgrid[2] - bgrid[1]   #step size on bgrid 
         k <- apply(p, 1, sum)*db     #vector of standardization constants 
          cdfba[,1] <<- db*p[,1]/k
         for ( j in 2:200)
             cdfba[,j] <<- cdfba[, j-1]+db*p[,j]/k
            
    
           
           NULL
     }
         
  ##############################################################
  #  Function rpost.82 
  #       
   #simulating a sample from the posterior using inverse cdf method
   # note:need to call init.rpost to setup cdfa & cdfba
   # see GCSR p 22/23 & p 84 for an explanation of the inverse cdf method 
  ##############################################################
   rpost.82 <- function(n) 
   {
    # Now generate p(a) and p(b |a) 
   u1<-runif(n)
   u2<- runif(n)    #uniform r.v.
   a<- rep(0,n)   #initialization
   b<-rep(0,n)     #initialization
   for(j in 1:n){
       k <- (1:200)[cdfa>=u1[j]][1]  # k= first index st cdfa[k]>u1 
       a[j] <-agrid[k]                         #a[j] =cdf-inv(u1[j]) 
       l<-(1:200)[cdfba[k,] >= u2[j]][1]     # l= first index st cdf[k]>u1 
       b[j]<- bgrid[l]                        #b[j]=cdfab-inv(u2[j])
       }
       return(cbind(a,b))  # return a (n,2) matrix of a,b vectors. 
       
     }
       
       #########################
       #  Main commands  
       #########################
       
       #marginal in alpha 
         p.agrid <<- apply(p,1,sum)  # note global assignment by by <<- 

       #Generate a posterior sample
      init.rpost.82()
      th <- rpost.82(200)    # simulate 200 draws from posterior 
      plot(th[,1],th[,2],xlab="ALPHA", 
ylab="BETA",xlim=c(-5,10),ylim=c(-10,40), pch=".")
      
 ###########################
 # Simulating LD95 
 ############################
 
       alpha = th[,1]
       beta = th[,2]
           
         lambda <- (log(0.95/0.05)-alpha)/beta 
        hist(lambda, col=0,xlab="LD95",ylab="P(LD95 | x)", prob=T )