多参数贝叶斯估计

logit<-function(x){
  y=log(x/(1-x))
  return(y)
}

bioassay<-data.frame(
x <- c(-0.86, -0.30, -0.05, 0.73),
n <- c(5, 5, 5, 5),
y <- c(0.01, 1, 3, 4.99),
r <- logit(y/n))
plot(x,r)
lm.bioassay<-lm(formula = r~x)
abline(lm.bioassay)
summary(lm.bioassay)

bioassay.post <- function(alpha = 0.1, beta = 5) {
  k <- 4
  x <- c(-0.86, -0.3, -0.05, 0.73)
  n <- c(5, 5, 5, 5)
  y <- c(0, 1, 3, 5)
  #prod连乘
  prod <- 1
  prod <- prod((exp(alpha + beta * x)/(1 + exp(alpha + beta * x)))^y *
                 (1/(1 + exp(alpha + beta * x)))^(n - y))
  return(prod)
}

mlpost <- function(alpha = 0.1, beta = 5) {
  -log(bioassay.post(alpha, beta))
}
library(bbmle)
mle(mlpost)


modedensity <- bioassay.post(0.87, 7.91)
alphax <- seq(-5, 10, length = 1000)
betay <- seq(-10, 40, length = 1000)
#Vectorize,能将一个不能进行向量化运算的函数进行转化，使之具备向量化运算功能。
post <- outer(alphax, betay, Vectorize(bioassay.post))
par(mfrow = c(1, 2))
contour(alphax, betay, post, levels = seq(0.05, 0.95, length = 10) * modedensity,
        xlim = c(-5, 10), ylim = c(-10, 40), xlab = quote(alpha), ylab = quote(beta),
        drawlabels = FALSE)

post <- post/sum(post)
posta <- apply(post, MARGIN = 1, FUN = sum)
w <- posta/sum(posta)
n <- 1000
ra <- rep(0, n)
rb <- rep(0, n)
for (j in 1:n) {
  ra[j] <- sample(alphax, 1, replace = T, prob = w)
  postb <- bioassay.post(ra[j], betay)
  wb <- postb/sum(postb)
  rb[j] <- sample(betay, 1, replace = T, prob = w)
}
plot(ra, rb, xlim = c(-5, 10), ylim = c(-10, 40), xlab = quote(alpha),
     ylab = quote(beta))

ld50 <- -ra/rb
hist(ld50,freq=FALSE,breaks=1000,xlim=c(-0.8,0.5),
       axes = TRUE, xlab="LD50",main="")