######################################################
#
#  Unit 2 - Example iv.
#
#  Sampling Importance Resampling
#
#  Normal data with known population and Cauchy 
#  prior for population mean.  
#
######################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
######################################################
# Sufficient statistics
# ---------------------
xbar  = 7.0
sig2n = 4.5

# Hyperparameters of the normal prior
# -----------------------------------
mu0   = 0.0
tau20 = 1.0

# Hyperparameters of the conjugate normal posterior
# -------------------------------------------------
tau21 = 1/(1/sig2n+1/tau20)
mu1   = tau21*(xbar/sig2n+mu0/tau20)

# Sampling from the prior
# -----------------------
set.seed(21386)
M  = 200
x1 = mu0+sqrt(tau20)*rt(M,1)

# Figure 3.2
# ----------
x          = seq(-14,14,length=5000)
prior      = dt((x-mu0)/sqrt(tau20),1)/sqrt(tau20)
likelihood = dnorm(x,xbar,sqrt(sig2n))
posterior  = prior*likelihood
w          = dnorm(x1,xbar,sqrt(sig2n))
w          = w/sum(w)

par(mfrow=c(2,1))
plot(x,posterior,xlab="",ylab="",main="",type="l",axes=F,ylim=c(0,0.004))
title("200 draws from the prior")
axis(1,at=c(-14,-5,0,5,14))
axis(2)
lines(x,likelihood/max(likelihood)*0.002,lty=2)
lines(x,prior/max(prior)*0.0036,lty=3)
for (i in 1:M){
  segments(x1[i],0,x1[i],-0.0001)
  segments(x1[i],0,x1[i],0.01*w[i],col=2)
}
abline(h=0)
text(0,0.0038,"PRIOR")
text(7,0.0022,"LIKELIHOOD")
text(4,0.0017,"POSTERIOR")

M   = 10000
x1  = mu0+sqrt(tau20)*rt(M,1)
x   = seq(-14,14,length=5000)
w   = dnorm(x1,xbar,sqrt(sig2n))
w   = w/sum(w)
ind = sample(1:M,size=M,replace=T,prob=w)
x2  = x1[ind]
hist(x2,xlab="",ylab="",prob=T,breaks=seq(min(x2),max(x2),length=20),xlim=c(-14,14),main="")
title("Posterior density \n based on 10000 draws")