###Poisson-gamma model
#y ~ Poisson(lambda)
#lambda ~ Gamma(alpha,beta) where E[lambda]=alpha*beta

#The integral that gives the marginal density f(y) can be solved in closed form,
#and establishes that f(y) is NegBin(alpha,1/(1+beta)) [see eqn (15.1)]


alpha=10; beta=1; y=15

###Evaluate f(y) using negative binomial
nb.prob=dnbinom(y,alpha,prob=1/(beta+1))
cat("\nProb that Y=",y,"is",nb.prob,"\n")

###Alternatively, do the integral numerical using the integrate function
#Define the integrand h(u) and log h(u)
log.const=-lgamma(y+1)-lgamma(alpha)-alpha*log(beta)
logh=function(u) -u+y*log(u)+(alpha-1)*log(u)-u/beta+log.const
h=function(u) exp(logh(u))
int.prob=integrate(h,lower=0.0001,upper=100)$value
cat("\nNumerical integration gives prob=",int.prob,"\n")

###Using the Laplace approximation
uhat=(y+alpha-1)*beta/(beta+1)
c.uhat=(y+alpha-1)/(uhat^2)
L.approx=h(uhat)*sqrt(2*pi/c.uhat)
cat("\nLaplace approximation gives prob =",L.approx,"\n")