# Example page 92 in Ed 2

# simulate values of lambda from the posterior distribution G(31/2, 6)
lambda = rgamma(10000, 31/2, rate=6)
hist(lambda, prob=T, xlab=expression(lambda), main="Posterior Distribution",
     ylim=c(0, max(dgamma(lambda, 31/2, rate=6))), breaks=50)
lines(density(lambda))
lines(sort(lambda), dgamma(sort(lambda), shape=31/2, rate=6), col=2, lty=2)

# need to sort the x values for the lines command above
# col = 2 is red
# lty = 2 is a dashed line

legend(3.5,.55, c("kernel density estimate", "exact Gamma density"), col=c(1,2), lty=c(1,2))

library(coda)
# use the function HPDinterval in the coda package to find an HPD interval
# coda assumes that output is in a format likel that produced by WinBUGS
# to use it when whe simulate our own posterior distribution we need to
# "coerce" the values into an MCMC object using the as.mcmc function

lambda = rgamma(100000, 31/2, rate=6)
HPDinterval(as.mcmc(lambda))

#> HPDinterval(as.mcmc(lambda))
#        lower    upper 
#var1 1.409210 3.940686
#attr(,"Probability")
#[1] 0.95
#
# compare to values from book
c(15.715, 45.452)/12
1.309583 3.787667

ypred = rpois(100000, lambda)