pdf(file="MCI-output.pdf",width=12,height=9)

###############################################################################
# Model: x|theta ~ N(theta,sigma^2)
###############################################################################
sigma = 40

###############################################################################
# Prior 1: theta ~ N(theta0,tau0^2)
###############################################################################
theta0 = 800
tau0   = 80

Mgrid      = 200
thetas     = seq(400,1200,length=Mgrid)
den.prior1 = dnorm(thetas,theta0,tau0)
plot(thetas,den.prior1,xlab="x or theta",ylab="Density",type="l",ylim=c(0,0.015),lwd=3)
title(paste("Model: x|theta ~ N(theta,",round(sigma^2,3),")\n Prior: theta ~ N(",theta0,",",round(tau0^2,3),")",sep=""))

# Predictive x ~ N(theta0,tau0^2+sigma^2)
mean  = theta0
stdev = sqrt(tau0^2+sigma^2)
xs     = seq(400,1200,length=Mgrid)
den.pred1 = dnorm(xs,mean,stdev)
lines(xs,den.pred1,col=2,lwd=3)

# Observated data
x = 975
points(x,0,pch=16,col=4,cex=2)
p1x = dnorm(x,mean,stdev)
tail1x = 1-pnorm(x,mean,stdev)

# Posterior
pi = sigma^2/(sigma^2+tau0^2)
stdev = sqrt(pi)*tau0
mean = pi*theta0+(1-pi)*x
den.post1 = dnorm(xs,mean,stdev)
lines(thetas,den.post1,col=6,lwd=2)

legend("topright",legend=c("p(theta)","p(x)",paste("p(theta|x=",x,")",sep="")),
col=c(1,2,6),lwd=3,bty="n")

legend("topleft",legend=c(
"Prior and likelihood conjugate, so all",
"derivations are done in closed form.",
"In general, approximate derivations are needed"))

legend("left",legend=c(
paste("p(x=",x,")=",round(p1x,6),sep=""),
paste("Pr(x>",x,")=",round(tail1x,6),sep="")))

###############################################################################
# Prior 2: theta ~ Student's t(theta0,tau0^2,nu0)
###############################################################################
nu0    = 4
tau0   = tau0*sqrt((nu0-2)/nu0)

den.prior2 = dt((thetas-theta0)/tau0,df=nu0)/tau0
plot(thetas,den.prior2,xlab="x or theta",ylab="Density",type="l",
     ylim=c(0,0.015),lwd=3)
title(paste("Model: x|theta ~ N(theta,",round(sigma^2,3),")\n Prior: theta ~ t(",theta0,",",round(tau0^2,3),",",nu0,")",sep=""))

# Predictive density (Monte Carlo integration)
set.seed(12345)
Mdraw      = 10000
theta.draw = theta0 + tau0*rt(M,df=nu0)
den.pred2 = rep(0,Mgrid)
for (i in 1:Mgrid)
  den.pred2[i] = mean(dnorm(xs[i],theta.draw,sigma))
lines(xs,den.pred2,col=2,lwd=3)
points(x,0,pch=16,col=4,cex=2)
x.draw = rnorm(Mdraw,theta.draw,sigma)
tail2x = mean(x.draw>x)

# Computing p(x=950) (via Monte Carlo integration)
p2x = mean(dnorm(x,theta.draw,sigma))

# Posterior
lines(thetas,den.prior2*dnorm(x,thetas,sigma)/p2x,col=6,lwd=2)
legend("topright",legend=c("p(theta)","p(x)",paste("p(theta|x=",x,")",sep="")),
col=c(1,2,6),lwd=3,bty="n")

legend("left",legend=c(
paste("p(x=",x,")=",round(p2x,6),sep=""),
paste("Pr(x>",x,")=",round(tail2x,6),sep="")))




###############################################################################
# Comparison
###############################################################################
plot(thetas,den.prior1,xlab="x or theta",ylab="Density",type="l",
     ylim=c(0,0.015),lwd=2)
lines(thetas,den.prior2,lwd=2,lty=2)
lines(xs,den.pred1,col=2,lwd=2)
lines(xs,den.pred2,col=2,lwd=2,lty=2)
lines(thetas,den.post1,col=6,lwd=2)
lines(thetas,den.prior2*dnorm(x,thetas,sigma)/p2x,col=6,lwd=2,lty=2)
legend("topright",legend=c("Gaussian prior","p(theta)","p(x)",paste("p(theta|x=",x,")",sep="")),
col=c(0,1,2,6),lwd=3,bty="n")
legend("topleft",legend=c("Student's t prior","p(theta)","p(x)",paste("p(theta|x=",x,")",sep="")),
col=c(0,1,2,6),lwd=3,lty=2,bty="n")
points(x,0,pch=16,col=4,cex=2)
legend("left",legend=c(
paste("p(x=",x,")=",round(p2x,6),sep=""),
paste("Pr(x>",x,")=",round(tail2x,6),sep="")))
legend("right",legend=c(
paste("p(x=",x,")=",round(p1x,6),sep=""),
paste("Pr(x>",x,")=",round(tail1x,6),sep="")))

BF = p1x/p2x
title(paste("Bayes factor=p(x|Gaussian)/p(x|Student's t)=",round(BF,3),sep=""))
dev.off()