logdt.hedi = function(x,m,s,nu){
  dt((x-m)/s,nu,log=TRUE)-log(s)
}
logdigamma.hedi = function(x,a,b){
  dgamma(1/x,a,b,log=TRUE)-2*log(x)
}
llike = function(y,theta,sig,nu){
  sum(logdt.hedi(y,theta,sig,nu))
}
lprior = function(theta,sig2,m0,sC0,a,b){
  dnorm(theta,m0,sC0,log=TRUE)+logdigamma.hedi(sig2,a,b)
}  
lpost = function(y,theta,sig2,nu,m0,sC0,a,b){
  sum(logdt.hedi(y,theta,sqrt(sig2),nu))+
  dnorm(theta,m0,sC0,log=TRUE)+
  logdigamma.hedi(sig2,a,b)
}
lprop = function(theta,sig2){
  dnorm(theta,4.5,3,log=TRUE)+dgamma(sig2,1.5,0.03,log=TRUE)
}


pdf(file="student-t.pdf",width=12,height=10)


#################################################################
# DATA SIMULATION 
#################################################################
set.seed(24435)
n     = 50
nu    = 4
theta = 10
sig   = sqrt(2.0)
y     = theta+sig*rt(n,df=nu)

breaks = seq(min(y),max(y),length=30)
hist(y,breaks=breaks,xlab="",main="",prob=TRUE)

#################################################################
# Back down to Earth where we only observe 
# y(1),...,y(n)
#################################################################

# Likelihood function and joint prior distribution
# theta ~ N(m0,C0) and sig2 ~ IG(n0/2,n0*s02/2)
m0  = 0
C0  = 0.5
n0  = 2
s02 = 2
sC0 = sqrt(C0)
a   = n0/2
b   = n0*s02/2

N        = 500
thetas   = seq(-2,12,length=N)
sig2s    = seq(0.001,70,length=N)
loglike  = matrix(0,N,N)
logprior = matrix(0,N,N) 
logpost  = matrix(0,N,N) 
logprop  = matrix(0,N,N) 
for (i in 1:N){
  for (j in 1:N){
    logprior[i,j] = lprior(thetas[i],sig2s[j],m0,sC0,a,b)    
    loglike[i,j]  = llike(y,thetas[i],sqrt(sig2s[j]),nu)
    logpost[i,j]  = lpost(y,thetas[i],sig2s[j],nu,m0,sC0,a,b)
    logprop[i,j]  = lprop(thetas[i],sig2s[j])
  }
}
prior = exp(logprior-max(logprior))
like  = exp(loglike-max(loglike))
post  = exp(logpost-max(logpost))
prop  = exp(logprop-max(logprop))

par(mfrow=c(1,2))
contour(thetas,sig2s,prior,col=4,drawlabels=FALSE,xlim=c(-2,2),
ylim=c(0,10),xlab=expression(theta),ylab=expression(sigma^2),main="Prior")
contour(thetas,sig2s,like,drawlabels=FALSE,xlim=c(8,12),ylim=c(0,5),
xlab=expression(theta),ylab=expression(sigma^2),main="Likelihood")

par(mfrow=c(1,1))
contour(thetas,sig2s,like,drawlabels=FALSE,xlab=expression(theta),
ylab=expression(sigma^2))
contour(thetas,sig2s,prior,col=4,add=TRUE,drawlabels=FALSE)
legend("topleft",legend=c("Likelihood","Prior"),col=c(1,4),lty=1)

par(mfrow=c(1,1))
contour(thetas,sig2s,like,drawlabels=FALSE,xlab=expression(theta),
ylab=expression(sigma^2))
contour(thetas,sig2s,prior,col=4,add=TRUE,drawlabels=FALSE)
contour(thetas,sig2s,post,col=6,add=TRUE,drawlabels=FALSE)
legend("topleft",legend=c("Likelihood","Prior","Posterior"),col=c(1,4,6),lty=1)

par(mfrow=c(1,1))
contour(thetas,sig2s,post,drawlabels=FALSE,xlab=expression(theta),
ylab=expression(sigma^2))
contour(thetas,sig2s,prop,col=2,add=TRUE)

############################################################
#  SAMPLING IMPORTANCE RESAMPLING (SIR)
############################################################
M = 20000
thetass = rnorm(M,4.5,3)
sig2ss  = rgamma(M,1.5,0.03)
w = rep(0,M)
for (i in 1:M)
  w[i] = lpost(y,thetass[i],sig2ss[i],nu,m0,sC0,a,b)-
         lprop(thetass[i],sig2ss[i])
w1 = exp(w-max(w))
ind = sample(1:M,size=M,replace=TRUE,prob=w1)
draws.sir = cbind(thetass[ind],sig2ss[ind])

par(mfrow=c(1,2))
plot(thetass,sig2ss,xlab=expression(theta),ylab=expression(sigma^2),col=grey(0.75),pch=16)
contour(thetas,sig2s,prop,col=2,add=TRUE,drawlabels=FALSE)
title("Proposal draws")
plot(draws.sir,xlab=expression(theta),ylab=expression(sigma^2),col=grey(0.75),pch=16)
contour(thetas,sig2s,post,col=2,add=TRUE,drawlabels=FALSE)
title("Posterior draws")

############################################################
# RANDOM WALK METROPOLIS-HASTINGS
############################################################
M0      = 100000
M       = 20000
LAG     = 1
niter   = M0+LAG*M

thetass2 = 5
sig2ss2  = 25
sd.t     = 0.1
sd.s2    = 1.0
draws.mh = matrix(0,niter,2)
set.seed(12345)
for (iter in 1:niter){
  draw.theta = rnorm(1,thetass2,sd.t)
  draw.sig2  = rnorm(1,sig2ss2,sd.s2)
  if(draw.sig2>0){
    logaccept = min(0,lpost(y,draw.theta,draw.sig2,nu,m0,sC0,a,b)-
                      lpost(y,thetass2,sig2ss2,nu,m0,sC0,a,b))
    if (log(runif(1))<logaccept){
      thetass2 = draw.theta
      sig2ss2  = draw.sig2
    }                
  }
  draws.mh[iter,1] = thetass2
  draws.mh[iter,2] = sig2ss2
}
draws.mh = draws.mh[seq(M0+1,niter,by=LAG),]

par(mfrow=c(2,3))
ts.plot(draws.mh[,1],xlab="Iterations",ylab="",main=expression(theta))
acf(draws.mh[,1],main="")
hist(draws.mh[,1],prob=TRUE,xlab="",main="")
ts.plot(draws.mh[,2],xlab="Iterations",ylab="",main=expression(sigma^2))
acf(draws.mh[,2],main="")
hist(draws.mh[,2],prob=TRUE,xlab="",main="")

par(mfrow=c(1,1))
contour(thetas,sig2s,prior,drawlabels=FALSE,xlab=expression(theta),
ylab=expression(sigma^2),col=4)
points(draws.mh,col=gray(0.8),pch=16,cex=0.5)
contour(thetas,sig2s,like,add=TRUE,drawlabels=FALSE)
contour(thetas,sig2s,post,col=6,add=TRUE,drawlabels=FALSE)
legend("topleft",legend=c("Likelihood","Prior","Posterior"),col=c(1,4,6),lty=1)

                  

############################################################
#  GIBBS SAMPLER WITH DATA AUGMENTATION
############################################################
n0s02   = n0*s02
n12     = (n0+n)/2
nu12    = (nu+1)/2

# initial values
sig2ss  = 30
thetass = 5

# running gibbs sampler
M0      = 100000
M       = 20000
LAG     = 1
niter   = M0+LAG*M
draws.gibbs = matrix(0,niter,2)
set.seed(12345)
for (iter in 1:niter){
  bs      = nu+(y-thetass)^2/sig2ss
  lambdas = 1/rgamma(n,nu12,bs/2)
  C1      = 1/(1/C0+sum(1/lambdas)/sig2ss)
  m1      = C1*(m0/C0+sum(y/lambdas)/sig2ss)
  thetass = rnorm(1,m1,sqrt(C1))
  n1s12   = n0s02+sum((y-thetass)^2/lambdas)
  sig2ss  = 1/rgamma(1,n12,n1s12/2) 
  draws.gibbs[iter,] = c(thetass,sig2ss)
}
draws.gibbs = draws.gibbs[seq(M0+1,niter,by=LAG),]

par(mfrow=c(2,3))
ts.plot(draws.gibbs[,1],xlab="Iterations",ylab="",main=expression(theta))
acf(draws.gibbs[,1],main="")
hist(draws.gibbs[,1],prob=TRUE,xlab="",main="")
ts.plot(draws.gibbs[,2],xlab="Iterations",ylab="",main=expression(sigma^2))
acf(draws.gibbs[,2],main="")
hist(draws.gibbs[,2],prob=TRUE,xlab="",main="")


par(mfrow=c(1,1))
contour(thetas,sig2s,prior,drawlabels=FALSE,xlab=expression(theta),
ylab=expression(sigma^2),col=4)
points(draws.gibbs,col=gray(0.8),pch=16,cex=0.5)
contour(thetas,sig2s,like,add=TRUE,drawlabels=FALSE)
contour(thetas,sig2s,post,col=6,add=TRUE,drawlabels=FALSE)
legend("topleft",legend=c("Likelihood","Prior","Posterior"),col=c(1,4,6),lty=1)


# Exact (grid approximation) versus SIR-based, RWMW-based and Gibbs-based approximations
h1 = thetas[2]-thetas[1]
mp.theta = apply(post,1,sum)
mp.theta = mp.theta/sum(mp.theta)/h1
h2 = sig2s[2]-sig2s[1]
mp.sig2 = apply(post,2,sum)
mp.sig2 = mp.sig2/sum(mp.sig2)/h2


par(mfrow=c(1,1))
plot(density(draws.gibbs[,1]),xlab="",main=expression(theta),lwd=2,lty=4)
lines(density(draws.sir[,1]),lwd=2,lty=2)
lines(density(draws.mh[,1]),lwd=2,lty=3)
lines(thetas,mp.theta,lwd=2,lty=1)
legend("topright",legend=c("Exact","SIR","RWMH","Gibbs"),lty=1:4,lwd=2)

par(mfrow=c(1,1))
plot(density(draws.gibbs[,2]),xlab="",main=expression(sigma^2),lwd=2,lty=4)
lines(density(draws.sir[,2]),lwd=2,lty=2)
lines(density(draws.mh[,2]),lwd=2,lty=3)
lines(sig2s,mp.sig2,lwd=2,lty=1)
legend("topright",legend=c("Exact","SIR","RWMH","Gibbs"),lty=1:4,lwd=2)

dev.off()

mtheta = sum(mp.theta*thetas*h)
mtheta1 = mean(draws.sir[,1])
mtheta2 = mean(draws.mh[,1])
mtheta3 = mean(draws.gibbs[,1])

neff.mh    = nrow(draws.mh)/(1+2*sum(acf(draws.mh[,1],plot=FALSE,lag=100)$acf))
neff.gibbs = nrow(draws.gibbs)/(1+2*sum(acf(draws.gibbs[,1],plot=FALSE,lag=100)$acf))
mcerror = c(0,sqrt(var(draws.sir[,1])/nrow(draws.sir)),
sqrt(var(draws.mh[,1])/neff.mh),
sqrt(var(draws.gibbs[,1])/neff.mh))

cbind(
c(mtheta,mtheta1,mtheta2,mtheta3),
c(0,nrow(draws.sir),nrow(draws.mh),nrow(draws.gibbs)),
round(c(0,nrow(draws.sir),neff.mh,neff.gibbs),0),
mcerror)