#Example:analysis of bioassay expriment

	
input <- function(LL=100){
  biossay <- data.frame(doses=c(-.863,-.296,-.053,.727),
	         rats=c(5,5,5,5),deaths=c(0.5,1,3,4.5),freq=c(0.5,1,3,4.5)/5)
  DD		   <- data.frame(y=log((biossay$freq)/(1-biossay$freq)),x=biossay$doses)
  estimates <- lm(y~x,data=DD)
  alpha.hat <- summary(estimates)$coeff[1,1]
  std.alpha <- summary(estimates)$coeff[1,2]
  beta.hat  <- summary(estimates)$coeff[2,1]
  std.beta  <- summary(estimates)$coeff[2,2]
  alpha       <- seq(-2,2,length=LL)
  beta       <-  seq(-5,10,length=LL)
  return(biossay,estimates,alpha,beta)
}

log.post<-function(alpha,beta,data=DD$biossay){
  ldens <- 0
  for (i in 1:length(data$doses)){
    theta <- 1/(1+exp(-alpha-beta*data$doses[i]))
    ldens  <- ldens + data$deaths[i]*log(theta)+(data$rats[i]-data$deaths[i])*log(1-theta)
  }
  ldens
}

plot.joint.post<-function(data){
  contours <- seq(.05,.95,.1)
  contour(DD$alpha,DD$beta,PP/max(PP),levels=contours,xlab="alpha",ylab="beta")
  mtext("Posterior density",3,line=1,cex=1.2)
}

grid.value<-function(data=DD$biossay,alpha,beta){
  ll    <- length(alpha)
  PP    <- matrix(NA,ll,ll)
  for(i in 1:ll){ 
    for(j in 1:ll){
      PP[i,j]<-exp(log.post(alpha[i],beta[j],data))	
    }}
  return(PP/sum(PP))
}


DD <- input(LL=200)
PP <- grid.value(DD$biossay,DD$alpha,DD$beta)
plot.joint.post(DD$bioassay)