logit <- function(p)  log(p/(1-p))

antilogit <- function(x) exp(x)/(1+exp(x))
#
# Binomial likelihood from a logistic regression model
likeli <- function(a,b,x,y,n) {
  p <- antilogit(a+b*x)
  exp(sum(y*log(p) + (n-y)*log(1-p)))
}
#
# Generate random values from a cdf on a grid
GenRandomFromGrid <- function(x, cdf.x) {
  u <- runif(1,0,1)
  indx <- u > cdf.x
  max(x[indx])
}

y <- c(0,1,3,5)
n <- rep(5,4)
x <- c(-0.86, -0.30, -0.05, 0.73)

EmpLogit <- log((y+0.5)/(n-y+0.5))
fit <- lm(EmpLogit ~ x)
summary(fit)

a.grid <- seq(from=-2, to=2, by=0.1)
b.grid <- seq(from=0, to=6, by=0.1)
n.a <- length(a.grid)
n.b <- length(b.grid)
z <- matrix(NA, nrow=n.a, ncol=n.b)
for (i in 1:n.a) {
  for (j in 1:n.b) {
    z[i,j] <- likeli(a.grid[i], b.grid[j], x, y, n)
  }
}

contour(a.grid, b.grid, z, xlab="alpha", ylab="beta", drawlabels=F)

a.grid <- seq(from=-4, to=10, by=0.1)
b.grid <- seq(from=-10, to=40, by=0.1)
n.a <- length(a.grid)
n.b <- length(b.grid)
z <- matrix(NA, nrow=n.a, ncol=n.b)
for (i in 1:n.a) {
  for (j in 1:n.b) {
    z[i,j] <- likeli(a.grid[i], b.grid[j], x, y, n)
  }
}
zz <- ifelse(is.nan(z), 0, z)  #  Get rid of NaNs from too small values; set them to 0.

contour(a.grid, b.grid, zz, xlab="alpha", ylab="beta", drawlabels=F)

p.ab <- zz/sum(zz)  # Normalize to get the joint posterior density on the grid.
contour(a.grid, b.grid, p.ab, xlab="alpha", ylab="beta", drawlabels=F)

p.a <- apply(p.ab, 1, sum)  # Get the marginal posterior density of alpha.
plot(a.grid, p.a/0.1, type="l", xlab="alpha", ylab="p(alpha | data)")

cdf.a <- cumsum(p.a)
plot(a.grid, cdf.a, type="l", xlab="x", ylab="Pr(alpha < x | data)")

#  Generate a random sample from the marginal posterior of alpha.
u <- runif(1000, 0, 1)  # Generate Unif[0,1] random variates.
a <- rep(NA, length(u))
for (i in 1:length(u)) {
  indx <- u[i] > cdf.a
  a[i] <- GenRandomFromGrid(a.grid, cdf.a)
}

hist(a, freq=F, xlab="alpha", ylim=c(0, 0.4), xlim=c(-4, 6), 
     breaks="Scott", main="Alpha")
lines(a.grid, p.a/0.1)

b <- rep(NA, length(a))
for (i in 1:length(a)) {
  indx <- (1:length(a.grid))[a.grid==a[i]]
  cdf.b.given.a <- cumsum(p.ab[indx,])/sum(p.ab[indx,])
  b[i] <- GenRandomFromGrid(b.grid, cdf.b.given.a)
}
plot(a,b,xlab="alpha", ylab="beta", pch=16, cex=0.5)

# If want to add some jitter
aa <- a + runif(length(a), -0.05, 0.05)
bb <- b + runif(length(b), -0.05, 0.05)
plot(aa,bb,xlab="alpha", ylab="beta", pch=16, cex=0.5)

LD50 <- -aa/bb
hist(LD50, xlab="LD50", freq=F, sub="(log dose scale)", breaks="Scott")