##############################################################################################################
# 
# Example 3.1 (pages 84) – Figure 3.1 (page 85) 
# 
# Poisson data (one single observation) with conjugate Gamma prior.
#
# Figure 3.1 exhibits exact posterior density (Gamma) as well as a normal 
#            approximation and two non-normal approximations.
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
alpha    = 3
beta     = 3
x        = 3
alpha1   = alpha + x
beta1    = beta  + 1
m        = (alpha1-1)/beta1
V        = (alpha1-1)/beta1^2
mphi     = log(alpha1/beta1)
Vphi     = 1/alpha1
mpsi     = sqrt((alpha1-1/2)/beta1)
Vpsi     = 1/(4*beta1)
xs       = seq(0.001,5,length=1000)
post     = matrix(0,1000,4)
post[,1] = dgamma(xs,alpha1,beta1)
post[,2] = dnorm(xs,m,sqrt(V))
post[,3] = dnorm(log(xs),mphi,sqrt(Vphi))/xs
post[,4] = dnorm(sqrt(xs),mpsi,sqrt(Vpsi))/(2*sqrt(xs))

# Figure 3.1
# ----------
par(mfrow=c(1,1))
plot(xs,post[,1],xlab="",ylab="",main="",axes=F,type="l",ylim=c(0,max(post)))
axis(1)
axis(2)
lines(xs,post[,2],lty=2)
lines(xs,post[,3],lty=3)
lines(xs,post[,4],lty=5)

# Monte Carlo approximation to the mean and variance of lambda when normal 
# approximations are derived for phi=log(lambda) and psi = sqrt(lambda) 
# ------------------------------------------------------------------------
xs   = rnorm(1000000,mphi,sqrt(Vphi))
xs   = exp(xs)
c(mean(xs),var(xs))
xs   = rnorm(1000000,mpsi,sqrt(Vpsi))
xs   = xs^2
c(mean(xs),var(xs))