## binomial model with non-conjugate prior
## ADM 6/15/2007
##
## all analysis is for y=6, n=10

## weird prior density
##
## input: p <- vector of abscissae
## output: fp <- density evaluated at each p
prior <- function(p) {
   fp <- matrix(NA, length(p),1)
   fp[(p <= 0.2 & p >= 0)] <- 0
   fp[(p <= 0.5 & p > 0.2)] <- 1/3
   fp[(p <= 0.9 & p > 0.5)] <- 2
   fp[(p <= 1   & p > 0.9)] <- 1
   return(fp)
}

## unnormalized posterior density
##
## inputs
##    p <- vector abscissae
##    y <- number of successes
##    n <- number of trials
##
## output: fpy <- unnormalized posterior density at each p
post.unnormalized <- function(p, y, n){
  fpy <- matrix(NA, length(p), 1)
  
  fpy[(p <= 0.2 & p >= 0)] <- 0

  indic <- (p <= 0.5 & p > 0.2) 
  fpy[indic] <- 1/3 * p[indic]^y * (1-p[indic])^(n-y)

  indic <- (p <= 0.9 & p > 0.5) 
  fpy[indic] <- 2 * p[indic]^y * (1-p[indic])^(n-y)

  indic <- (p <= 1 & p > 0.9)
  fpy[indic] <- 1  * p[indic]^y * (1-p[indic])^(n-y)
  
  return(fpy)
}

## plot prior (must do these piecewise so lines are not connected)
ruler1 <- 0:200/1000
ruler2 <- 201:500/1000
ruler3 <- 501:900/1000
ruler4 <- 901:1000/1000
plot(ruler1, prior(ruler1), type="l", lwd=4, col="grey40", xlim=c(0,1), ylim=c(0,3),
     xlab="pi", ylab="f(pi)")
lines(ruler2, prior(ruler2), lwd=4, col="grey40")
lines(ruler3, prior(ruler3), lwd=4, col="grey40")
lines(ruler4, prior(ruler4), lwd=4, col="grey40")

## plot unnormalized posterior
plot(ruler1, post.unnormalized(ruler1, 6, 10), type="l", lwd=4, col="black", 
   xlim=c(0,1), ylim=c(0,0.0025), xlab="pi", ylab="f(pi|y)")
lines(ruler2, post.unnormalized(ruler2, 6, 10), lwd=4, col="black")
lines(ruler3, post.unnormalized(ruler3, 6, 10), lwd=4, col="black")
lines(ruler4, post.unnormalized(ruler4, 6, 10), lwd=4, col="black")

## find normalizing constant
cat("\n## Normalizing Constant\n")
constant <- integrate(post.unnormalized, 0, 1, y=6, n=10)
print(constant)

## plot prior and normalized posterior on the same graph
plot(ruler1, prior(ruler1), type="l", lwd=4, col="grey40", xlim=c(0,1), 
   ylim=c(0,4), xlab="pi", ylab="")
lines(ruler2, prior(ruler2), lwd=4, col="grey40")
lines(ruler3, prior(ruler3), lwd=4, col="grey40")
lines(ruler4, prior(ruler4), lwd=4, col="grey40")
lines(ruler1, post.unnormalized(ruler1, 6, 10)/constant[[1]], lwd=4, 
   col="black")
lines(ruler2, post.unnormalized(ruler2, 6, 10)/constant[[1]], lwd=4, 
   col="black")
lines(ruler3, post.unnormalized(ruler3, 6, 10)/constant[[1]], lwd=4, 
   col="black")
lines(ruler4, post.unnormalized(ruler4, 6, 10)/constant[[1]], lwd=4, 
   col="black")
legend(0, 4, legend=c("Prior", "Posterior"),
       col=c("grey40", "black"), lwd=4)

## draw from the posterior
ruler <- seq(0, 1, 0.001)
probs <- post.unnormalized(ruler, 6, 10)
probs <- probs/sum(probs)
draws <- sample(ruler, 10000, prob=probs, replace=TRUE)
cat("\n## Posterior Mean\n")
print(mean(draws))
cat("\n## Posterior Quantiles\n")
print(quantile(draws))
cat("\n## Pr(pi > 0.4)\n")
print(sum(draws > 0.4)/10000)
cat("\n")