#############################################################################################
#
# Posterior inference via 
#
#  1) Sampling Importance Resampling (SIR)
#  2) Metropolis-Hastings (MCMC)
#
# Source of the banana-shape joint posterior:
#
# Recursive Pathways to Marginal Likelihood Estimation with Prior-Sensitivity Analysis
# By Ewan Cameron and Anthony Pettitt
# Statistical Science
# 2014, Vol. 29, No. 3, 397–419
# DOI: 10.1214/13-STS465
# Institute of Mathematical Statistics, 2014
#
#############################################################################################
library(MASS)
install.packages("shape")
library(shape)


##########################################
# Bivariate posterior
# Prior of theta1 is U(-0.5,1)
# Prior of theta2 is U(-0.5,1)
# Prior independence of theta1 and theta2
##########################################
prior.likelihood = function(theta1,theta2){
  prior=(1/1.5)^2
  likelihood=exp(-(10*(0.45-theta1))^2/4-(20*(theta2/2-theta1^4))^2)
  return(prior*likelihood)
}

######################## 
# SIR
######################## 

time.sir = system.time({
# Step 1:  Sampling from proposal
M = 100000
Mr = M/10
theta1.prior = runif(M,-0.5,1)
theta2.prior = runif(M,-0.5,1)

# Step 2: Computing resampling weights
w = rep(0,M)
for (i in 1:M)
  w[i] = prior.likelihood(theta1.prior[i],theta2.prior[i])
  
# Step 3: Resampling
ind = sample(1:M,size=Mr,replace=TRUE,prob=w)
theta1.post = theta1.prior[ind]
theta2.post = theta2.prior[ind]
thetas.sir = cbind(theta1.post,theta2.post)
})


par(mfrow=c(2,2))
plot(theta1.prior,theta2.prior,xlab=expression(theta[1]),ylab=expression(theta[2]),
     main="Prior draws",pch=".")
plot(theta1.post,theta2.post,xlim=range(theta1.prior),ylim=range(theta2.prior),
     xlab=expression(theta[1]),ylab=expression(theta[2]),pch=".",
     main="Posterior draws")
hist(theta1.post,prob=TRUE,xlab=expression(theta[1]),main="",xlim=c(-0.5,1),col=3)
title("Marginal posterior")
abline(h=1/1.5,col=2,lwd=3)
hist(theta2.post,prob=TRUE,xlab=expression(theta[2]),main="",xlim=c(-0.5,1),col=3)
title("Marginal posterior")
abline(h=1/1.5,col=2,lwd=3)


# Posterior means
E1.sir = mean(theta1.post)
E2.sir = mean(theta2.post)


######################## 
# Metropolis-Hastings
######################## 


# Initial values
theta1 = 0.25
theta2 = 0.25
M0     = 10000
M      = 10000
L      = 9
niter  = M0+M*L
thetas.mh = matrix(0,niter,2) 
time.mh = system.time({
for (i in 1:niter){
  theta1d = runif(1,-0.5,1)
  theta2d = runif(1,-0.5,1)
  alpha = min(1,prior.likelihood(theta1d,theta2d)/prior.likelihood(theta1,theta2))
  if (runif(1)<alpha){
    theta1 = theta1d
    theta2 = theta2d
  }
  thetas.mh[i,] = c(theta1,theta2)
}
thetas.mh = thetas.mh[seq(M0+1,niter,by=L),]
})

par(mfrow=c(2,3))
ts.plot(thetas.mh[,1],xlab="Iterations",ylab="draws",main=expression(theta[1]))
acf(thetas.mh[,1],main="")
plot(density(thetas.mh[,1]),xlab="",main="")
lines(density(thetas.sir[,1]),col=2)
legend("topright",legend=c("SIR","MCMC"),col=2:1,lty=1,bty="n")
ts.plot(thetas.mh[,2],xlab="Iterations",ylab="draws",main=expression(theta[2]))
acf(thetas.mh[,2],main="")
plot(density(thetas.mh[,2]),xlab="",main="")
lines(density(thetas.sir[,2]),col=2)


z.sir =  kde2d(thetas.sir[,1], thetas.sir[,2], n=100)
z.mh  =  kde2d(thetas.mh[,1], thetas.mh[,2], n=100)
par(mfrow=c(1,2))
plot(thetas.sir,cex=0.5,xlab=expression(theta[1]),ylab=expression(theta[2]))
points(thetas.mh,col=2,cex=0.5)
legend("topleft",legend=c("SIR","MCMC"),col=1:2,pch=16)
contour(z.sir, drawlabels=FALSE,nlevels=20,xlab=expression(theta[1]),ylab=expression(theta[2]))
contour(z.mh, drawlabels=FALSE,add=TRUE,col=2,nlevels=20)



E1.mh = mean(thetas.mh[,1])
E2.mh = mean(thetas.mh[,2])

c(E1.sir,E1.mh)
c(E2.sir,E2.mh)

time.sir
time.mh