############################################################################################
#
# Example 3.7 (pages 102-3) – Figure 3.4 (page 104)  - Weighted resampling algorithm (SIR)
#
############################################################################################
#
# 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)}


# 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
# --------------------------------------
M   =  100
th1 =  seq(-20,50,length=M)
th2 =  seq(-2,6,length=M)
f   =  matrix(0,M,M)
for (i in 1:M)
  for (j in 1:M)
    f[i,j] = post(c(th1[i],th2[j]))
f = f/sum(f)

# Sampling importance resampling (SIR) algorithm
# ----------------------------------------------
set.seed(9023749)
M       = 10000
ths     = matrix(0,M,2)
ths[,1] = -20+70*runif(M)
ths[,2] = -2+8*runif(M)
w = NULL
for (i in 1:M)
  w = c(w,post(c(ths[i,1],ths[i,2])))
w    = w/sum(w)
ind  = sample(1:M,size=M,replace=T,prob=w)
ths1 = ths[ind,]

# Posterior summaries
# -------------------
matrix(c(apply(ths1,2,mean),
sqrt(apply(ths1,2,var)),
apply(ths1,2,quantile025),
apply(ths1,2,quantile975)),2,4)

# Curve fitting
# -------------
N  = 20
x1 = seq(min(x),max(x),length=N)
f1 = matrix(0,M,N)
for (i in 1:N)
  f1[,i] = ths1[,1]*(1-exp(-x1[i]*ths1[,2]))
q050 = apply(f1,2,median)
q025 = apply(f1,2,quantile025)
q975 = apply(f1,2,quantile975)

# Figure 3.4
# ----------
par(mfrow=c(1,2))
plot(ths1,xlab=expression(theta[1]),ylab=expression(theta[2]),pch=".",xlim=range(th1),ylim=range(th2),axes=F)
axis(1)
axis(2)
title("(a)")
contour(th1,th2,f,drawlabels=F,nlevels=7,axes=F,levels=c(0.00005,0.0005,seq(0.001,0.005,by=0.001)),add=T)
for (i in 1:M)
   points(ths1[i,1],ths1[i,2],col=1,pch=".")
plot(x,y,ylim=range(f1),xlab="time (days)",ylab="BOD (mg/liter)",main="",axes=F)
axis(1)
axis(2)
title("(b)")
lines(x1,q025,lty=2)
lines(x1,q050)
lines(x1,q975,lty=2)