# Normalized posterior for model 1
post1 = function(x1,x2,mu){
  dnorm(mu,10/21*(1+x1+x2),sqrt(10/21))	
}

# Unnormalized posterior for model 2
post2 = function(x1,x2,mu){
  dnorm(mu,0,sqrt(10))/((1+(mu-x1)^2)*(1+(mu-x2)^2))
}

# Observed data
x = c(-4.3,3.2)

# Ploting the posterior densities
mu = seq(-10,10,length=1000)
par(mfrow=c(1,2))
plot(mu,post1(x[1],x[2],mu),type="l",lwd=2,ylab="Prior*likelihood")
title("Model 1: Gaussian")
points(x[1],0,pch=16,col=3)
points(x[2],0,pch=16,col=3)

plot(mu,post2(x[1],x[2],mu),type="l",lwd=2,ylab="Prior*likelihood")
title("Model 2: Cauchy")
points(x[1],0,pch=16,col=3)
points(x[2],0,pch=16,col=3)

# Computing via simple Reimann integration the 
# normalizing constant for model M2
M      = 100000
mus   = seq(-10,10,length=M)
h       = mus[2]-mus[1]
prod = rep(0,M)
for (i in 1:M)
prod[i] = h*post2(x[1],x[2],mus[i])
kappa = sum(prod)


par(mfrow=c(1,1))
plot(mu,post1(x[1],x[2],mu),type="l",lwd=2,xlab=expression(mu),
ylab="Posterior density")
lines(mu,post2(x[1],x[2],mu)/kappa,col=2,lwd=2)
lines(mu,dnorm(mu,0,sqrt(10)),col=4,lwd=2)
legend("topleft",legend=c("Gaussian prior for mu","Gaussian model","Cauchy model"),col=c(4,1:2),lwd=2,lty=1)
points(x[1],0,pch=16,col=3)
points(x[2],0,pch=16,col=3)




# Using SIR to transform draws from the prior into draws from the posterior
# for model M2 (Cauchy model)
N = 1000000
mus = rnorm(N,0,sqrt(10))
w = 1/((1+(mus-x[1])^2)*(1+(mus-x[2])^2))
ind = sample(1:N,size=N,prob=w,replace=TRUE)
mus1 = mus[ind]
hist(mus1,prob=TRUE,xlab=expression(mu),ylab="Posterior density",main="",
ylim=c(0,0.2))
lines(mu,post2(x[1],x[2],mu)/kappa,col=2,lwd=2)
legend("topleft",legend=c("Exact posterior","MC-based posterior draws"),col=2:1,lty=1,lwd=2)