##########################################
#### Quantile based Credible Intervals ###
##########################################


### Example 1, Chapter 2 notes: (example 2.1, page 46-47 in book)

# Data for the Cabinet Duration example:

N <- c(38,28,27,20,17,15,15,15,15,14,12)
# The sum of the vector N gives the total number of observations
dur <- c(0.833,1.070,1.234,1.671,2.065,2.080,2.114,2.168,2.274,2.629,2.637)
# The sum of the vector N*dur gives the sum of all the durations (of all the x_i's)

# A 90% quantile-based credible interval for theta (reciprocal of mean duration):
qgamma(0.05, shape=sum(N), rate=sum(N*dur))
qgamma(0.95, shape=sum(N), rate=sum(N*dur))

### Example 2, Chapter 2 notes:

x <- 2; n <- 10

# A 95% quantile-based credible interval for theta:

qbeta( c(.025, .975), x+1, n-x+1)


######################
#### HPD Intervals ###
######################

#### Example 1, Chapter 2 notes: (example 2.2, page 51 in book)

N <- c(38,28,27,20,17,15,15,15,15,14,12)
# The sum of the vector N gives the total number of observations
dur <- c(0.833,1.070,1.234,1.671,2.065,2.080,2.114,2.168,2.274,2.629,2.637)
# The sum of the vector N*dur gives the sum of all the durations (of all the x_i's)

ruler <- seq(0.45, 0.75, length=5000)
# The "ruler" is basically the horizontal line segment that will be "dropped" onto the posterior density
# The "0.45, 0.75" indicates where the segment should begin and end

g.vals <- round(dgamma(ruler, shape=sum(N), rate=sum(N*dur)), 2)
start <- 1000; stop <- length(g.vals)
target <- 0.90; tol <- 0.005
# The "target" is the desired probability level for the interval.

done <- FALSE; i <- start
while (i < stop & done == FALSE) {
 j <- length(g.vals)/2
 while (j <= stop & done == FALSE) {
  if (g.vals[i] == g.vals[j])  {
   L <- pgamma(ruler[i], shape=sum(N), rate=sum(N*dur))
   H <- pgamma(ruler[j], shape=sum(N), rate=sum(N*dur))
   if (((H-L)<(target+tol)) & ((H-L)>(target-tol)))
    done <- TRUE
  }
  j <- j+1
 }
      i <- i+1
 
}
k <- g.vals[i]; HPD.L <- ruler[i]; HPD.U <- ruler[j]

print(paste(target*100, "% HPD interval:", HPD.L, "to", HPD.U))

### Example 2, Chapter 2 notes:

library(TeachingDemos)
# May need to install the TeachingDemos package first?
# If so, type at the command line:  install.packages("TeachingDemos", dependencies=F)
# while plugged in to the internet.

x <- 2; n <- 10

hpd(qbeta, shape1 = x+1, shape2= n-x+1, conf=0.95)

# Using hpd function for example 1 above:

hpd(qgamma, shape=sum(N), rate=sum(N*dur), conf=0.90)

###################################
#### More With Conjugate Priors ###
###################################

# Anthropology Example (example 2.8, page 70-71 of book)

par(oma=c(1,1,1,1), mar=c(0,0,0,0), mfrow=c(2,1))

x <- c(1,1,1,1,0,1,1,0,1,0,1,1,1,0,1,1,1,1,1,1,0,0,0,1)
y<- sum(x)
n <- length(x)

my.seq <- seq(0,1,length=300)

########### Analysis with beta(15,2) prior:

a<-15; b<-2

beta.prior <- dbeta(my.seq, a, b) 
# Defining the prior
beta.posterior <- dbeta(my.seq, y+a, n-y+b)
# Defining the posterior

plot(my.seq, beta.prior, ylim=c(-0.7,9.5), yaxt='n', 
     xlab="Values of p", ylab="Probability distribution", type='l', lty=3)
# Plotting the prior
lines(my.seq, beta.posterior)
# Plotting the posterior

### Point estimates for p:

posterior.mean <- (y+a)/(a+b+n)

posterior.median <- qbeta(0.50, y+a, n-y+b)

posterior.mode <- my.seq[beta.posterior==max(beta.posterior)]

print(paste("posterior.mean=", round(posterior.mean,3), 
      "posterior.median=", round(posterior.median,3), "posterior.mode=", round(posterior.mode,3)))

### Interval estimate for p:

hpd.95 <- hpd(qbeta, shape1=y+a, shape2=n-y+b, conf=0.95)
segments(hpd.95[1], 0, hpd.95[2], 0, lwd=4)
# Showing the 95% HPD interval on the graph

text(mean(hpd.95), -0.4, "95% HPD Interval", cex=0.6)
text(0.25,5,paste("beta(", a, ",", b, ") prior:", 
     "95% HPD Interval at: [", round(hpd.95[1], 3), ",", round(hpd.95[2],3), "]"), cex=0.8)


########### Analysis with beta(1,1) (i.e., Uniform) prior:

a<-1; b<-1

beta.prior <- dbeta(my.seq, a, b) 
# Defining the prior
beta.posterior <- dbeta(my.seq, y+a, n-y+b)
# Defining the posterior

plot(my.seq, beta.prior, ylim=c(-0.7,9.5), yaxt='n', 
     xlab="Values of p", ylab="Probability distribution", type='l', lty=3)
# Plotting the prior
lines(my.seq, beta.posterior)
# Plotting the posterior

### Point estimates for p:

posterior.mean <- (y+a)/(a+b+n)

posterior.median <- qbeta(0.50, y+a, n-y+b)

posterior.mode <- my.seq[beta.posterior==max(beta.posterior)]

print(paste("posterior.mean=", round(posterior.mean,3), 
      "posterior.median=", round(posterior.median,3), "posterior.mode=", round(posterior.mode,3)))

### Interval estimate for p:

hpd.95 <- hpd(qbeta, shape1=y+a, shape2=n-y+b, conf=0.95)
segments(hpd.95[1], 0, hpd.95[2], 0, lwd=4)
# Showing the 95% HPD interval on the graph

text(mean(hpd.95), -0.4, "95% HPD Interval", cex=0.6)
text(0.25,5,paste("beta(", a, ",", b, ") prior:", 
     "95% HPD Interval at: [", round(hpd.95[1], 3), ",", round(hpd.95[2],3), "]"), cex=0.8)


par(mfrow=c(1,1)) # resetting plotting window