#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter5: Section5.2.3
#
# Poisson models with conjugate priors
#
# Author : Tomohiro Ando
#
################################################################################


# Data generation

x <- c(rep(0,len=6),1,4)

n <- 8


# Setting hyperparameters

#alpha <- beta <- 2 #Model 1

alpha <- beta <- 10  #Model 2


# Posterior inference

alphahat <- alpha+sum(x)

betahat  <- beta+n


# Prior density, Likelihood, Posterior density

l.range <- seq(0,4,len=100)

Prior <- dgamma(l.range, shape=alpha, scale=1/beta)

Posterior <- dgamma(l.range, shape=alphahat, scale=1/betahat)

Likelihood <-  exp( sum(x)*log(l.range)-n*l.range+log(24) )


# Marginal likelihood calculation

Marginallikelihood <- gamma(alphahat)/gamma(alpha)*beta^alpha/(betahat)^alphahat/24

print(Marginallikelihood)



#Plot

plot(l.range,Prior,type="l",xlab="",ylab="",bty="L", ylim=c(0,2.5))
lines(l.range,Posterior,lty=2)

par(new=TRUE)
plot(l.range, Likelihood, type="l",lty=1,axes=FALSE, xlab="", ylab="")