####################################################################################################
#
#
#  Normal model and normal prior - univariate example
#  
#  Goal: illustrate the potential fallacy of using the prior as proposal in the SIR algorithm
#  
#  By Hedibert Freitas Lopes
#  This version: May 22nd 2021
#
####################################################################################################

# Gaussian data and Gaussian prior
# x1,...,xn|theta ~ N(theta,1)
# theta ~ N(0,10)
n = 5
x = c(3,3,4,5,5)

# Closed form (analytical/exact) posterior inference
# theta|x1,...xn ~ N(theta1,tau12)
xbar   = mean(x)
theta1 = (10/(1/n+10))*xbar
tau12 = 10*(1/n)/(1/n+10)
curve(dnorm(x,0,sqrt(10)),from=-10,to=10,n=1000,col=2,lwd=2,
      ylab="Density",xlab=expression(theta),ylim=c(0,1))
curve(dnorm(x,theta1,sqrt(tau12)),from=-3,to=7,n=1000,lwd=2,add=TRUE)
curve(dnorm(x,xbar,sqrt(1/n)),from=-3,to=7,n=1000,col=4,lwd=2,add=TRUE)
legend("topleft",legend=c("Prior","Likelihood","Posterior"),col=c(2,4,1),
        bty="n",lty=1,lwd=2)

# The following posterior summaries 
# can easily be computed analytically
# A. Pr(theta>4.5|x1,...,x5) = c
# B. E(theta|x1,...,x5) = a
# C. Pr(theta<L|x1,...,x5) = 0.025
# D. Pr(theta<U|x1,...,x5) = 0.975

A = 1-pnorm(4.5,theta1,sqrt(tau12))
B = theta1
C = qnorm(0.025,theta1,sqrt(tau12))
D = qnorm(0.975,theta1,sqrt(tau12))

summary = c(A,B,C,D)


# What if we were not capable of computing A,B,C and D exactly?
# -------------------------------------------------------------
likelihood = function(theta){
  exp(-0.5*(n*theta^2-2*theta*n*xbar))
}
prior = function(theta){
  dnorm(theta,0,sqrt(10))
}

M = 10000
N = 10000

# SIR: Prior as proposal
theta.draw1 = rnorm(M,0,sqrt(10))
w           = likelihood(theta.draw1)*prior(theta.draw1)/prior(theta.draw1)
theta.post1 = sample(theta.draw1,size=N,replace=TRUE,prob=w)

# SIR: Uniform proposal in the range (0,10)
theta.draw2 = runif(M,0,10)
w           = likelihood(theta.draw2)*prior(theta.draw2)/(1/10)
theta.post2 = sample(theta.draw2,size=N,replace=TRUE,prob=w)

par(mfrow=c(2,2))
hist(theta.draw1,prob=TRUE,xlab="theta",main="Proposal draws\nProposal equal to prior",xlim=c(-20,20))
hist(theta.post1,prob=TRUE,xlab="theta",main="Posterior draws",xlim=c(0,10))
hist(theta.draw2,prob=TRUE,xlab="theta",main="Proposal draws\nProposal NOT equal to prior",xlim=c(-20,20))
hist(theta.post2,prob=TRUE,xlab="theta",main="Posterior draws",xlim=c(0,10))


summary1 = c(mean(theta.post1>4.5),
mean(theta.post1),
quantile(theta.post1,0.025),
quantile(theta.post1,0.975))


summary2 = c(mean(theta.post2>4.5),
mean(theta.post2),
quantile(theta.post2,0.025),
quantile(theta.post2,0.975))

cbind(summary,summary1,summary2)







# What is the size of the MC approximation error?
Rep      = 100
summary1 = matrix(0,Rep,4)
summary2 = matrix(0,Rep,4)
for (rep in 1:Rep){
theta.draw1 = rnorm(M,0,sqrt(10))
w           = likelihood(theta.draw1)*prior(theta.draw1)/prior(theta.draw1)
theta.post1 = sample(theta.draw1,size=N,replace=TRUE,prob=w)
theta.draw2 = runif(M,0,10)
w           = likelihood(theta.draw2)*prior(theta.draw2)/(1/10)
theta.post2 = sample(theta.draw2,size=N,replace=TRUE,prob=w)
summary1[rep,] = c(mean(theta.post1>4.5),
mean(theta.post1),
quantile(theta.post1,0.025),
quantile(theta.post1,0.975))
summary2[rep,] = c(mean(theta.post2>4.5),
mean(theta.post2),
quantile(theta.post2,0.025),
quantile(theta.post2,0.975))
}



names = c("A. Pr(theta>4.5|x1,...,x5)","B. E(theta|x1,...,x5)",
"C. Pr(theta<L|x1,...,x5) = 0.025","D. Pr(theta<U|x1,...,x5) = 0.975")
par(mfrow=c(2,2))
for (i in 1:4){
  boxplot(summary1[,i],summary2[,i],names=c("Proposal 1","Proposal 2"))
  title(names[i])
  abline(h=summary[i],col=2)
}