#############################################################################################################
#
# Example 6.4 Part 2: (pages 216-7) – Figures 6.9 (page 218)  - Centered regressor
#
# In this example time to failure (y) of motorettes were tested at different temperatures.  
# 
# A simple linear regression is fit
#
#         t[i] ~ N(alpha + beta z[i],sigma2)
#
#  z[i] = 1000/(x[i]+273.2)
#  t[i] = log10(y[i])
#
# sigma2 = 0.2592
# 
# Prior: p(alpha,beta) = c
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
loglike = function(m,n,y,x,a,b,s){
  max(sum(dnorm(y[1:m],a+b*x[1:m],s,log=T))+
  sum(log(1-pnorm(y[(m+1):n],a+b*x[(m+1):n],s))),-5000)
}

y = c(1764,2772,3444,3542,3780,4860,5196,408,408,1344,1344,1440,408,408,504,504,504,
      rep(8064,10),rep(5448,3),rep(1680,5),rep(528,5))
x = c(rep(170,7),rep(190,5),rep(220,5),rep(150,10),rep(170,3),rep(190,5),rep(220,5))
y = log10(y)
x = 1000/(x+273.2)
mx = mean(x)
x = x - mx
n = 40
m = 17

s  = 0.2592
N  = 100
as = seq(3.2,3.8,length=N)
bs = seq(2.8,6.2,length=N)
like = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
    like[i,j] = loglike(m,n,y,x,as[i],bs[j],s)


par(mfrow=c(1,2))
plot(x,y,xlab="Temperatures (Celsius)",ylab="Time to failures (hours)")
contour(as,bs,exp(like),xlab="alpha",ylab="beta",drawlabels=F)

################################################################################################################################
#                                                    SIMULATIONS                                                               #
################################################################################################################################
a0     = c(3.3,3.7,3.7)
b0     = c(3,3,5.5)
sd     = 0.1
nn     = length(a0)
niter0 = 200
niter  = 20000
SIMUL  = array(0,c(3,3,niter+1,2))

# [1] Random Walk Metropolis (single move)
# ----------------------------------------
set.seed(9237497)
draw  = array(0,c(3,niter+1,2))
for (j in 1:nn){
  a  = a0[j]
  b  = b0[j]
  draws=c(a,b)
  for (iter in 1:niter){
    a1 = rnorm(1,a,sd)
    accept = min(0,loglike(m,n,y,x,a1,b,s)-loglike(m,n,y,x,a,b,s))
    if (log(runif(1))<accept){
      a=a1
    }
    b1 = rnorm(1,b,sd)
    accept = min(0,loglike(m,n,y,x,a,b1,s)-loglike(m,n,y,x,a,b,s))
    if (log(runif(1))<accept){
      b=b1
    }
    draws = rbind(draws,c(a,b))
  }
  draw[j,,] = draws
}
SIMUL[1,,,] = draw

# [2] Random Walk Metropolis (block move)
# ---------------------------------------
set.seed(9237497)
draw  = array(0,c(3,niter+1,2))
for (j in 1:3){
  a  = a0[j]
  b  = b0[j]
  draws=c(a,b)
  for (iter in 1:niter){
    a1 = rnorm(1,a,sd)
    b1 = rnorm(1,b,sd)
    accept = min(0,loglike(m,n,y,x,a1,b1,s)-loglike(m,n,y,x,a,b,s))
    if (log(runif(1))<accept){
      a=a1
      b=b1
    }
    draws = rbind(draws,c(a,b))
  }
  draw[j,,] = draws
}
SIMUL[2,,,] = draw

# [3] Independence Metropolis algorithm
# -------------------------------------
reg   = lm(y[1:m]~x[1:m])
X     = cbind(1,x[1:m])
mean  = reg$coef
s     = 0.2592
V     = s*solve(t(X)%*%X)
cV    = t(chol(V))
icV   = solve(cV)

set.seed(9237497)
draw  = array(0,c(3,niter+1,2))
for (j in 1:3){
  a  = a0[j]
  b  = b0[j]
  z  = icV%*%(matrix(c(a,b),2,1)-mean)
  draws=c(a,b)
  for (iter in 1:niter){
    z1     = rnorm(2)
    draw1  = mean + cV%*%z1
    a1     = draw1[1]
    b1     = draw1[2]
    accept = min(0,loglike(m,n,y,x,a1,b1,s)+sum(dnorm(z,log=T))-loglike(m,n,y,x,a,b,s)-sum(dnorm(z1,log=T)))
    if (log(runif(1))<accept){
      a = a1
      b = b1
      z = z1
    }
    draws = rbind(draws,c(a,b))
  }
  draw[j,,] = draws
}
SIMUL[3,,,] = draw


# Comparing all algorithms - Figure 6.9
# -------------------------------------
met = c("Random walk Metropolis (single)","Random walk Metropolis (block)","Independence Metropolis")
par(mfrow=c(3,3))
for (k in 1:3){
  draw = SIMUL[k,,,]
  if (k<4){
    contour(as,bs,exp(like),xlab=expression(beta[0]),ylab=expression(beta[1]),drawlabels=F,xlim=range(as),ylim=range(bs))
    title(met[k])
    for (j in 1:3){
      points(a0[j],b0[j],pch=16)
      for (i in 1:min(niter0,niter)){
        segments(draw[j,i,1],draw[j,i,2],draw[j,i+1,1],draw[j,i,2])
        segments(draw[j,i+1,1],draw[j,i,2],draw[j,i+1,1],draw[j,i+1,2])
      }
    }
  }else{
    contour(as,bs,exp(like),xlab=expression(beta[0]),ylab=expression(beta[1]),drawlabels=F,xlim=c(3.3,3.7),ylim=c(3,5.5))
    title(met[k])
    for (j in 1:3)
      points(draw[j,,])
  }
  ts.plot(draw[1,,1],xlab="iterations",ylab="",main=expression(beta[0]),ylim=range(as))
  for (i in 2:3) 
    lines(draw[i,,1])
  ts.plot(draw[1,,2],xlab="iterations",ylab="",main=expression(beta[1]),ylim=range(bs))
  for (i in 2:3) 
    lines(draw[i,,2])
}

# Posterior summary
# -----------------
stats = NULL
for (i in 1:3)
  for (j in 1:2)
    stats = rbind(stats,c(mean(SIMUL[i,,,j]),median(SIMUL[i,,,j]),quantile(SIMUL[i,,,j],0.025),quantile(SIMUL[i,,,j],0.975)))

# Effective sample sizes (Table 6.3)
# ----------------------------------
M    = 3*(niter + 1) 
neff = matrix(0,3,2)
for (i in 1:3)
  for (k in 1:2){
    series = c(SIMUL[i,1,,k],SIMUL[i,2,,k],SIMUL[i,3,,k])
    neff[i,k] = M/(1+0.01*sum((200-2:200)*acf(series,lag=200,plot=F)$acf[2:200]))
  }