###############################################
 #Binomial Data Example
 #Placenta Previa (Example 2.5 in GRCS)
 #1/4/06
 ###############################################
 
 ####
 ##Summarizing the Exact Posterior Distribution of theta -- Beta(438,544)
 ####
 
 post.mean <- 0.446                                 #posterior mean
 post.sd <- 0.016                                   #posterior sd
 
 #Posterior Quantiles
 post.median <- qbeta(.5,438,544)                   #median
 post.interval.95 <- qbeta(c(0.025,0.975),438,544)  #95% posterior interval
 
 #Plot of the Posterior Density Function
 x <- seq(0,1,length=500)
 plot(x,dbeta(x,438,544),type='l',xlab="theta",ylab="p(theta|y)",ylim=c(0,25))
 abline(v=post.mean,col="red")
 abline(v=jitter(post.median,.1),col="blue")    #jitter the median so that it is visible
 abline(v=post.interval.95,col="green")
 legend(.6,10,c("mean","median","95% interval"),col=c("red","blue","green"),text.col="black",lty = 1)
 
 #QUESTION:  How much evidence is there that the proportion of female placenta previa births is less than 0.485?
 plot(x,dbeta(x,438,544),type='l',xlab="theta",ylab="p(theta|y)",ylim=c(0,25),xlim=c(.3,.7))
 abline(v=0.485,col="red")
 
 pbeta(0.485,438,544)  #area less that 0.485

 ####
 ##Approximating the Posterior Distribution of theta
 ####
 
 #Get 1000 draws from the Beta(438,544) distribution
 theta <-  rbeta(1000,438,544)
 
 #Sample summaries
 post.mean.approx <- mean(theta)      #approximate posterior mean
 post.sd.approx <- sd(theta)          #approximate posterior sd
 post.median.approx <- median(theta)  #approximate posterior median
 post.interval.95.approx <- sort(theta)[c(25,975)]  #approximate 95% posterior interval
 
 #Histogram of the samples
 hist(theta,prob=T,ylim=c(0,25))
 lines(x,dbeta(x,438,544),col="red")
 
 #Compare approximation using different sample sizes
 par(mfrow=c(3,1))    # 3 plots in one window
 num.draws <- c(100,1000,10000)
 plot.titles <- c("Histogram -- 100 Draws","Histogram -- 1000 Draws","Histogram -- 10,000 Draws")
 for(i in 1:3){
    theta <-  rbeta(num.draws[i],438,544)
    hist(theta,prob=T,ylim=c(0,30),main=plot.titles[i])
    lines(x,dbeta(x,438,544),col="red")
    }
    
 ####
 ##Normal approximation to the posterior
 ####

 par(mfrow=c(1,1))
 x <- seq(0,1,length=500)
 plot(x,dbeta(x,438,544),type='l',xlab="theta",ylab="p(theta|y)",ylim=c(0,25),xlim=c(.3,.7))
 lines(x,dnorm(x,post.mean,post.sd),col='red')
 
 ####
 ##Transformations of the posterior samples
 ####
 
 par(mfrow=c(2,2))
 theta <-  rbeta(1000,438,544)
 hist(theta,main="Histogram -- theta",prob=T)
 hist(log(theta),main="Histogram -- log(theta)",prob=T)
 hist(theta/(1-theta),main="Histogram of the sex ratio",prob=T)
 hist(log(theta/(1-theta)),main="Histogram -- logit(theta)",prob=T)