#################################################################################
#
#  Unit 5 - Example i.
#
#     Several methods to compute the normalizing constant f(y)
#       
#################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
#################################################################################
prior  = function(theta){dt((theta-theta0)/sqrt(tau02),1)/tau02}
like   = function(theta){dnorm(theta,xbar,sqrt(sig2n))}
post   = function(theta){prior(theta)*like(theta)}
g      = function(theta){dnorm(theta,th.hat,sqrt(V))}
alpha  = function(theta){1/sqrt(g(theta)*post(theta))}
g0     = function(theta){prior(theta)}
gk     = function(theta){post(theta)}
gi     = function(theta,i){g0(theta)^(1-lambda[i])*gk(theta)^(lambda[i])}

# Situation of example 2.1 and 2.3 of Gamerman and Lopes (2006)
# -------------------------------------------------------------
theta0 = 0.0
tau02  = 1.0
xbar   = 7.0
sig2n  = 4.5
theta  = seq(-15,15,length=1000)

# Riemman's integration
# ---------------------
thetas = seq(-15,15,length=10000000)
f.true = (thetas[2]-thetas[1])*sum(post(thetas))


# Finding the mode of the posterior density
# -----------------------------------------
func   = function(theta){-log(post(theta))}
th.hat = nlm(func,0.05)$estimate

# Sampling from the posterior using a SIR algorithm
# -------------------------------------------------
set.seed(192477)
M   = 100000
th  = rnorm(M,th.hat,2)
w   = post(th)/dnorm(th,th.hat,2)
w   = w/sum(w)
ind = sample(1:M,size=M,replace=T,prob=w)
th1 = th[ind]
V   = var(th1)   # Variance of the normal approximation

# f0: Normal approximation estimator
# ----------------------------------
f0 = post(th.hat)*sqrt(2*pi)*sqrt(V)

# f1: simple Monte Carlo estimator
# --------------------------------
set.seed(192477)
th = theta0+sqrt(tau02)*rt(M,1)
f1 = mean(like(th))

# f4: harmonic mean estimator
# ---------------------------
f4 = 1/mean(1/like(th1))

# f5: Newton and Raftery
# ----------------------
set.seed(192477)
deltas = seq(0.01,0.1,by=0.01)
th     = theta0+sqrt(tau02)*rt(M,1)
u      = runif(M)
f5s    = NULL
for (delta in deltas){
  M1  = sum(u<delta)
  th2 = c(th[1:M1],th1[1:(M-M1)])

  f5 = f4 
  for (i in 1:10){
    num = sum(like(th2)/(delta*f5+(1-delta)*like(th2)))
    den = sum(1/(delta*f5+(1-delta)*like(th2)))
    f5  = num/den
  }
  f5s = c(f5s,f5)
}

# f6: Gelfand and Dey estimator
# -----------------------------
f6 = 1/mean(g(th1)/post(th1))

# f7: Annealed importance sampling
# --------------------------------
set.seed(192477)
M      = 5000
K      = 20
sd     = 1
lambda = c(seq(0.1,0.5,length=K/4),seq(0.51,1.0,length=3*K/4))
thetas = NULL
ws     = NULL
for (j in 1:M){
  th     = theta0
  ths    = th
  thstar = th+sqrt(tau02)*rt(1,1)
  alpha  = min(1,gi(thstar,1)/gi(th,1)*g0(th)/g0(thstar))
  th     = ifelse(runif(1)<alpha,thstar,th)
  w      = gi(th,1)/g0(th)
  ths    = c(ths,th)
  for (i in 2:K){
    thstar = rnorm(1,th,sd)
    alpha  = min(1,gi(thstar,i)/gi(th,i))
    th     = ifelse(runif(1)<alpha,thstar,th)
    w      = w*gi(th,i)/gi(th,i-1)
    ths    = c(ths,th)
  }
  thetas = rbind(thetas,ths)
  ws = c(ws,w)
}
f7 = sum(ws)/M

# f8: optimal bridge estimator (with geometric bridge estimator as initial value)
# -------------------------------------------------------------------------------
alpha = function(theta){1/sqrt(g(theta)*post(theta))}
set.seed(192477)
th3 = rnorm(M,th.hat,sqrt(V))
fbs = mean(alpha(th3)*post(th3))/mean(alpha(th1)*g(th1))
w   = post(th1)/g(th1)
wt  = post(th3)/g(th3)
s1  = length(th1)/(length(th1)+length(th2))
s2  = length(th2)/(length(th1)+length(th2))
f8  = fbs
for (i in 1:10)
  f8 = mean(wt/(s1*wt+s2*f8))/mean(1/(s1*w+s2*f8))

# f11: Chib and Jeliazkov estimator
# ---------------------------------
set.seed(192477)
f11s = NULL
phis = c(4.80,4.85,4.90,4.95,5.00,5.05,5.076856,5.10,5.15,
         5.20,5.25,5.30,5.35,5.384173,5.40,5.45,5.50,seq(5.55,6.0,by=0.05))
for (phi in phis){
  th4  = rnorm(M,phi,sqrt(V))
  accept1 = post(phi)/post(th1)*dnorm(th1,phi,sqrt(V))/dnorm(phi,th1,sqrt(V))
  accept1[accept1>1]=1
  accept2  = post(th4)/post(phi)*dnorm(phi,th4,sqrt(V))/dnorm(th4,phi,sqrt(V))
  accept1[accept2>1]=1
  pphi = mean(accept1*dnorm(phi,th1,sqrt(V)))/mean(accept2)
  f11  = like(phi)*prior(phi)/pphi
  f11s = c(f11s,f11)
}

# Figure 7.1  (from Gamerman and Lopes, 2006): prior, likelihood, posterior and normal approximation
# ---------------------------------------------------------------------------------------------------
par(mfrow=c(2,2))
breaks = seq(min(th1),max(th1),length=50)
plot(theta,prior(theta),xlab="",ylab="",main="(a)",ylim=c(0,0.32),type="l",axes=F)
axis(1)
axis(2)
lines(theta,like(theta),lty=2)
lines(theta,dnorm(theta,th.hat,sqrt(V)),lty=3)
hist(th1,xlab="",ylab="",main="(b)",breaks=breaks,prob=T)
lines(theta,post(theta)/f.true)
plot(deltas,10000*f5s,type="b",xlab=expression(delta),ylab="10000 f(y)",main="(c)",axes=F,ylim=10000*range(f5s,f11s))
axis(1,at=deltas)
axis(2)
abline(h=10000*f.true)
plot(phis,10000*f11s,type="b",xlab=expression(phi),ylab="10000 f(y)",main="(d)",axes=F,ylim=10000*range(f5s,f11s))
axis(1)
axis(2)
abline(h=10000*f.true)
points(phis[7],f11s[7],pch=16)
points(phis[14],f11s[14],pch=15)

# Table 7.2  (from Gamerman and Lopes, 2006)
# ------------------------------------------
fff   = round(c(f0,f1,f4,f5s[10],f6,f7,f8,f11s[14]),8)
error = round(100*abs(fff-f.true)/f.true,2)
cbind(fff,error)