##############################################################################################################
#
# Example 7.11 (pages 277) – Figure 7.6 (page 278)
#
# Multiple try Metropolis and Random walk Metropolis algorithms to sample from 
#
# This example is related to Example 3.7 (pages 102-3).
#
##############################################################################################################
#
# NONLINEAR MODEL MODEL WITH EXPONENTIAL DECAY WITH A FIXED RATE CONSTANT
#
# Biochemical oxygen demand (BOD)
#
# Tanner (1996, page 169)
# Marske (1967) 
# Bates and Watts (1988)
#
# x = time (days)
# y = BOD (mg/Liter)
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
quantile025 = function(x){quantile(x,0.025)}
quantile975 = function(x){quantile(x,0.975)}
post = function(th){sum((y-th[1]*(1-exp(-x*th[2])))^2)^(-2)*dnorm(th[1],20,20)*dnorm(th[2],0,1.5)}

# Data
# ----
x = c(1,2,3,4,5,7)
y = c(8.3,10.3,19,16,15.6,19.8)

# Contours of the posterior distribution
# --------------------------------------
th1 =  seq(0,60,length=100)
th2 =  seq(0,7,length=100)
f   =  matrix(0,100,100)
for (i in 1:100)
  for (j in 1:100)
    f[i,j] = post(c(th1[i],th2[j]))
f = f/sum(f)


# Random Walk Metropolis 
# ----------------------
set.seed(9023749)
M       = 10000
ths     = matrix(0,M+1,2)
sd      = c(25,3)
th      = c(10,0)
ths[1,] = th
for (iter in 1:M){
    # sampling theta1 
    theta = rnorm(1,th[1],sd[1])
    accept = min(1,post(c(theta,th[2]))/post(th))
    if (runif(1)<accept){
      th[1] = theta
    }
    # sampling theta2 
    theta = rnorm(1,th[2],sd[2])
    accept = min(1,post(c(th[1],theta))/post(th))
    if (runif(1)<accept){
      th[2] = theta
    }
    ths[iter+1,] = th
}

# Random Walk Multiple Try Metropolis 
# -----------------------------------
k        = 10
M1       = 1000
ths1     = matrix(0,M1+1,2)
th1      = matrix(0,k,2)
sd       = c(25,3)
th       = c(10,0)
ths1[1,] = th
w       = rep(0,k)
for (iter in 1:M1){
    # sampling theta1 
    theta = rnorm(k,th[1],sd[1])
    for (i in 1:k)
      w[i] = post(c(theta[i],th[2]))
    w   = w/sum(w)
    th2 = sample(theta,size=1,prob=w)
    th3 = c(rnorm(k-1,th2,sd[1]),th2)
    num = 0.0
    den = 0.0
    for (i in 1:k){
      num = num + post(c(theta[i],th[2]))
      den = den + post(c(th3[i],th[2]))
    }
    accept = min(1,num/den)
    if (runif(1)<accept){
      th[1] = th2
    }
    # sampling theta2 
    theta = rnorm(k,th[2],sd[2])
    for (i in 1:k)
      w[i] = post(c(th[1],theta[i]))
    w   = w/sum(w)
    th2 = sample(theta,size=1,prob=w)
    th3 = c(rnorm(k-1,th2,sd[2]),th2)
    num = 0.0
    den = 0.0
    for (i in 1:k){
      num = num + post(c(th[1],theta[i]))
      den = den + post(c(th[1],th3[i]))
    }
    accept = min(1,num/den)
    if (runif(1)<accept){
      th[2] = th2
    }
    ths1[iter+1,] = th
}

# Figure 7.6
# ----------
th1 = seq(0,60,length=100)
th2 = seq(0,7,length=100)
th1 = seq(-1,70,length=100)
th2 = seq(-0.12,7,length=100)

par(mfrow=c(2,3))

plot(ths[,1],xlab="iterations",ylab=expression(theta[1]),main="(a)",ylim=range(th1),axes=F,type="l")
axis(1);axis(2)
plot(ths[,2],xlab="iterations",ylab=expression(theta[2]),main="(b)",ylim=range(th2),axes=F,type="l")
axis(1);axis(2)
plot(ths,xlab=expression(theta[1]),ylab=expression(theta[2]),pch=".",xlim=range(th1),ylim=range(th2),axes=F)
title("(c)")
axis(1)
axis(2)
contour(th1,th2,f,nlevels=7,levels=c(0.0002,0.0004,0.0008,0.0015,0.003,0.008),axes=F,drawlabels=F,add=T)
for (i in 1:M)
   points(ths[i,1],ths[i,2],col=1,pch=".",cex=3)

plot(ths1[,1],xlab="iterations",ylab=expression(theta[1]),main="(d)",ylim=range(th1),axes=F,type="l")
axis(1);axis(2)
plot(ths1[,2],xlab="iterations",ylab=expression(theta[2]),main="(e)",ylim=range(th2),axes=F,type="l")
axis(1);axis(2)
plot(ths1,xlab=expression(theta[1]),ylab=expression(theta[2]),pch=".",xlim=range(th1),ylim=range(th2),axes=F)
title("(f)")
axis(1);axis(2)
contour(th1,th2,f,nlevels=7,levels=c(0.0002,0.0004,0.0008,0.0015,0.003,0.008),axes=F,drawlabels=F,add=T)
for (i in 1:M1)
   points(ths1[i,1],ths1[i,2],col=1,pch=".",cex=3)

# Statistical summary
# -------------------
round(c(mean(ths[,1]),sqrt(var(ths[,1])),quantile(ths[,1],0.025),quantile(ths[,1],0.975)),3)
round(c(mean(ths[,2]),sqrt(var(ths[,2])),quantile(ths[,2],0.025),quantile(ths[,2],0.975)),3)
round(c(mean(ths1[,1]),sqrt(var(ths1[,1])),quantile(ths1[,1],0.025),quantile(ths1[,1],0.975)),3)
round(c(mean(ths1[,2]),sqrt(var(ths1[,2])),quantile(ths1[,2],0.025),quantile(ths1[,2],0.975)),3)

c(M/(1+2*sum(acf(ths[,1],lag=200,plot=F)$acf[2:200])),
M/(1+2*sum(acf(ths1[,1],lag=200,plot=F)$acf[2:200])),
M1/(1+2*sum(acf(ths[,2],lag=200,plot=F)$acf[2:200])),
M1/(1+2*sum(acf(ths1[,2],lag=200,plot=F)$acf[2:200])))

c(length(unique(ths[,1]))/M,length(unique(ths[,2]))/M,length(unique(ths1[,1]))/M1,length(unique(ths1[,2]))/M1)