## code to accompany the binomial example(s) from lecture 3
## of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy adapted from Peter Hoff

## plotting the posterior
theta <- seq(0,1,length=1000)
plot(theta, dbeta(theta, 119, 12), type="l", lwd=2,
     ylab="density", xlab="theta", bty="n")
abline(h=1, col="gray", lwd=2)
legend(x="topleft", bty="n", c("posterior", "prior"), lwd=2,
       col=c("black", "gray"))

## add confidence region to happiness
q <- qbeta(c(0.025, 0.975), 119, 12)
abline(v=q, , lwd=2, lty=2, col=2)
legend("left", bty="n", c("95% CI"), lwd=2, lty=2, col=2)

## simple confidence region for a beta distribution
a <- 1; b <- 1  ## prior
n <- 10; y <- 2 ## data
q <- qbeta(c(0.025, 0.975), a+y, b+n-y)
q

## plot the confidence region
plot(theta, dbeta(theta, a+y, b+n-y), type="l", lwd=2,
     ylab="density", xlab="theta", bty="n")
abline(v=q, lwd=2, lty=2, col=2)
legend(x="topright", bty="n", c("posterior", "95% CI"),
       lwd=2, lty=1:2, col=1:2)


## finding the hpd region (numerically) via Hoff
db <- dbeta(theta, a+y, b+n-y)
ord <- order(-db)
xpx <- cbind(theta[ord], db[ord])
xpx <- cbind(xpx, cumsum(xpx[,2])/sum(xpx[,2]))

## hpd:
##
## extract the HPD at prob p from the first
## two columns of xpx, which are theta and p(theta),
## as defined above

hpd <- function(x, dx, p)
  {
    ## normalize the probability density,
    ## and put in decreasing order
    px <- dx/sum(dx)
    pxs <- sort(px, decreasing=TRUE)

    ## find the smallest point with CDF less than p
    ct <- min(pxs[cumsum(pxs) < p])

    ## done, return the range of thetas greater than ct
    return(range(x[px>=ct]))
  }

## initialize HPD plot
plot(theta, dbeta(theta, a+y, b+n-y), type="l", lwd=2,
     ylab="density", xlab="theta", bty="n")

## add legend
legend("topright", c("95% q-CI", "50% HPD","75% HPD","95% HPD"), 
       col=2:5, lty=2:5, lwd=2, bty="n")

## add 95% quantile lines
lines(x=c(q[1], q[1], q[2], q[2]),
      y=dbeta(c(0, q[1], q[2], 0) ,a+y, b+n-y),
      col=2, lwd=2 ,lty=2)

## add 50% HPD lines
hpd50 <- hpd(xpx[,1],xpx[,2],0.5)
lines(x=c(hpd50[1], hpd50[1],  hpd50[2], hpd50[2]),
      y=dbeta(c(0,hpd50[1], hpd50[2],0), a+y, b+n-y),
      col=3, lty=3, lwd=2)

## add 75% HPD lines
hpd75 <- hpd(xpx[,1], xpx[,2], 0.75)
lines(x=c(hpd75[1], hpd75[1], hpd75[2], hpd75[2]),
      y=dbeta(c(0, hpd75[1], hpd75[2], 0), a+y, b+n-y),
      col=4, lty=4, lwd=2)

## add 95% HPD lines
hpd95 <- hpd(xpx[,1],xpx[,2],0.95)
lines(x=c(hpd95[1], hpd95[1], hpd95[2], hpd95[2]),
      y=dbeta(c(0, hpd95[1], hpd95[2], 0), a+y, b+n-y),
      col=5, lty=5, lwd=2)
hpd95


## the Laplace Parisian births example
pbeta(0.5, 241945+1, 251527+1, lower.tail=FALSE)

## placenta previa example
pbeta(0.485, 437+1, 543+1, lower.tail=FALSE)

## plotting the full posterior 
plot(theta, dbeta(theta, 438, 544), type="l", lwd=2,
     ylab="density", xlab="theta", bty="n")
q <- qbeta(c(0.025, 0.975), 438, 544)
abline(v=q, lwd=2, lty=2, col=2)
legend(x="topright", bty="n", c("posterior", "95% CI"),
       lwd=2, lty=1:2, col=1:2)
points(x=0.45, y=0, col="pink", cex=2)