#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter5: Section5.3.1
#
# Bernoulli distribution with conjugate prior
#
# Author : Tomohiro Ando
#
################################################################################

# Data generation

x <- c(rep(0,len=8),rep(1,len=2))

n <- 10


# Setting hyperparameters

#alpha <- 0.1 #Model 1
#alpha <- 2 #Model 2
#alpha <- 4 #Model 3
alpha <- 8 #Model 4

beta  <- 4


# Posterior inference

alphahat <- alpha+sum(x)

betahat  <- beta+n-sum(x)


# Prior density, Likelihood, Posterior density

p.range <- seq(0.01,0.99,len=100)

Prior <- dbeta(p.range, shape1=alpha, shape2=beta)

Posterior <- dbeta(p.range, shape1=alphahat, shape2=betahat)

Likelihood <- dbinom(sum(x),size=n,prob=p.range) 


# Marginal likelihood calculation

Marginallikelihood <- choose(10,2)*gamma(alpha+beta)/(gamma(alpha)*gamma(beta))*gamma(alphahat)*gamma(betahat)/(gamma(n+alpha+beta))

print(Marginallikelihood)


#Plot

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

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