#Fitting a logistic model for a bioassay experiment
setwd("~/opetus/bayesian_analysis/luennot")
library(rjags)

x <- c(-0.86,-0.30,-0.05,0.73)
y <- c(0,1,3,5)
n <- c(5,5,5,5)
xpred <- seq(-1,1,by=0.1)
data <- list(k=4,n=n,x=x,y=y,xpred=xpred,K=len(xpred))
bioassay.jag <- jags.model("Bioassay.txt",data,inits=list(beta=1))
bioassay.coda <- coda.samples(bioassay.jag,c("alpha","beta","LD50","theta.pred"),10000)
a <- summary(bioassay.coda)

#Quantiles for the prediction
med <- a$quantiles[-(1:3),3]
low <- a$quantiles[-(1:3),1]
upp <- a$quantiles[-(1:3),5]

#postscript("fig12.ps", horizontal = FALSE, onefile = FALSE,
#          paper = "special", width = 10, height = 6)
op <- par(mfrow=c(1,2))
plot(xpred,upp,type="l",lty=3,ylim=c(0,1),xlab="dose",ylab="probability of death")
lines(xpred,low,lty=3)
lines(xpred,med)
points(x,y/n,pch=19)
hist(bioassay.coda[[1]][,"LD50"],main="",xlab="LD50",breaks=40)

#dev.off()