#############################################################################################
#
# Posterior inference via Sampling Importance Resampling (SIR)
#
# 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
#
#############################################################################################


####################################
# One-dimensional problem
# Prior of theta is U(-0.5,1)
####################################
prior.likelihood = function(theta){
  prior=1/1.5
  likelihood=exp(-(10*(0.45-theta))^2/4-(20*(0.065-theta^4))^2)
  return(prior*likelihood)
}

N     = 10000
theta = seq(-0.5,1,length=N)
post  = prior.likelihood(theta)
h     = theta[2]-theta[1]
py    = sum(post*h)
post  = post/py

# Expectation, variance and standard deviation
E = sum(theta*post*h)
E2 = sum(theta^2*post*h) 
V = E2 - E^2
S = sqrt(V)
c(E,S)

plot(theta,post,type="l",xlab=expression(theta),main="Posterior density",ylab="Density")
abline(v=E,col=2)


# SIR
# Step 1:  Sampling from proposal
M      = 1000
thetas = runif(M,-0.5,1)

# Step 2: Computing resampling weights
w      = prior.likelihood(thetas)
w      = w/sum(w)
plot(thetas,w,pch=16,cex=0.5,xlab=expression(theta),ylab="Weights")

# Step 3: Resampling
thetas.post = sample(thetas,size=M,replace=TRUE,prob=w)


E.sir = mean(thetas.post)
S.sir = sqrt(var(thetas.post))

c(E,E.sir)
c(S,S.sir)

hist(thetas.post,prob=TRUE,breaks=seq(-0.5,1,length=40),xlab=expression(theta),main="")
lines(theta,post,col=2,lwd=2)
title("Posterior: grid vs SIR")


# How good (or bad) is the Monte Carlo approximation to the posterior mean?
Ms = c(seq(10,100,by=10),seq(200,1000,by=1000),
       seq(2000,10000,by=1000),seq(20000,100000,by=10000),seq(200000,1000000,by=100000))
nM = length(Ms)
E.sir = rep(0,nM)
for (i in 1:nM){
  thetas = runif(Ms[i],-0.5,1)
  w      = prior.likelihood(thetas)
  thetas.post = sample(thetas,size=M,replace=TRUE,prob=w)
  E.sir[i] = mean(thetas.post)
}


plot(log10(Ms),E.sir,ylab="Posterior mean",xlab="Monte Carlo size (log10)",type="b",pch=16)
abline(h=E,col=2,lwd=2)
title("Monte Carlo approximation of posterior mean")





##########################################
# 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)
}

N      = 200
theta1 = seq(-0.5,1,length=N)
theta2 = seq(-0.5,1,length=N)

post = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
    post[i,j] = prior.likelihood(theta1[i],theta2[j])

h = theta1[2]-theta1[1]
py = sum(apply(post*h^2,1,sum))
post = post/py

E1 = sum(t(post)%*%theta1*h^2)
E2 = sum(post%*%theta2*h^2)

c(E1,E2,py)

par(mfrow=c(1,1))
contour(theta1,theta2,post,nlevels=15,xlab=expression(theta[1]),ylab=expression(theta[2]),
drawlabels=FALSE,main="Posterior")


# SIR

# Step 1:  Sampling from proposal
M = 100000
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=M,replace=TRUE,prob=w)
theta1.post = theta1.prior[ind]
theta2.post = theta2.prior[ind]

# Computing normalizing constant (prior predictive)
py.mc = mean(w)


par(mfrow=c(2,2))
plot(theta1.prior,theta2.prior,xlab=expression(theta[1]),ylab=expression(theta[2]),
     main="Prior draws vs posterior evaluations",pch=".")
contour(theta1,theta2,post,nlevels=15,drawlabels=FALSE,add=TRUE,col=3,lwd=2)
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 vs posterior evaluations")
contour(theta1,theta2,post,nlevels=15,drawlabels=FALSE,add=TRUE,col=3,lwd=2)
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)



# Marginal priors and posteriors
tab = rbind(
summary(theta1.prior),
summary(theta2.prior),
summary(theta1.post),
summary(theta2.post))
rownames(tab) = c("Prior theta1","Prior theta2","Posterior theta1","Posterior theta2")
round(tab,3)


E1.mc = mean(theta1.post)
E2.mc = mean(theta2.post)


c(E1,E1.mc)
c(E2,E2.mc)
c(py,py.mc)