# Count data
set.seed(12356)
n  = 10
lambda.true = 2
y  = rpois(n,lambda.true)
sn = sum(y)

mean(y)

var(y)

plot(table(y)/n)

###########################################################################################
# Model 0: 
# IID Poisson data with rate parameter lambda: yi|lambda ~ Poi(lambda)
# Gamma prior for lambda: lambda|a0,b0 ~ Gamma(a0,b0)
###########################################################################################

a0 = 4
b0 = 2

# Prior predictive density
# ------------------------
term1 = log(1/prod(factorial(y)))
term2 = log(b0^a0/gamma(a0))
term3 = lgamma(a0+sn)-(a0+sn)*log(b0+n)
logpred0 = term1+term2+term2
pred0=exp(logpred0)
c(logpred0,pred0)

# Posterior density
# -----------------
a1 = a0+sn
b1 = b0+n


lambda     = seq(0.001,10,length=1000)
plot(lambda,dgamma(lambda,a1,b1),col=1,type="l",lwd=2,ylab="Density")
lines(lambda,dgamma(lambda,a0,b0),col=2,lwd=2)
legend("topright",legend=c("Prior","Posterior"),col=2:1,lwd=2)
points(sn/n,0,col=3,pch=16,cex=2)

###########################################################################################
# Model 1: 
# IID Poisson data with rate parameter lambda: yi|lambda ~ Poi(lambda)
# Log-normal prior for lambda: lambda|mu0,sig0 ~ Log-normal(mu0,sig0^2)
###########################################################################################
mu0  = 0.4904
sig0 = sqrt(0.405465)

term1 = -sum(log(factorial(y)))
term2 = -0.5*log(2*pi*sig0^2)

# Riemann-type integral approximation 
h     = 0.005
theta = seq(1,2.5,by=h)
term3    = log(sum(theta^(sn-1)*exp(-n*theta)*exp(-0.5/sig0^2*(log(theta)-mu0)^2)*h))

logpred1 = term1+term2+term3
pred1    = exp(logpred1)
c(logpred1,pred1)



###########################################################################################
# Comparison
###########################################################################################
fgamma     = dgamma(lambda,a0,b0)
flognormal = dnorm(log(lambda),mu0,sig0)/lambda

plot(lambda,fgamma,xlab="lambda",ylab="Density",type="l",lty=1,lwd=3,ylim=c(0,0.5))
lines(lambda,flognormal,col=2,lwd=3)
legend("topright",legend=c("Gamma prior","Log-normal prior"),col=1:2,lwd=3)
points(sn/n,0,col=3,pch=16,cex=2)

c(pred0,pred1)
c(logpred0,logpred1)

pred0/pred1


######################################################
# How about trying a MC Integration approach?
######################################################