################################################################################
#
# Third homework assignment - Solution
# PhD in Business Economics 
# Advanced Bayesian Econometrics
# Hedibert Freitas Lopes 
# Due at 9am, February 18th, 2021.
#
# Nonlinear Gaussian regression
# 
# Let us consider the context of Example 6.1 (pages 192-194) of 
# Gamerman and Lopes (2006), where the response variable y, is 
# the velocity of an enzymatic reaction (in counts/min/min) and 
# the regressor, x, is substrate concentration (in ppm), where 
# the enzyme has been treated with Puromycin.  Check the book 
# webpage at http://www.dme.ufrj.br/mcmc/chapter6.html. 
# 
# Reference: Carlin and Louis (1998) Bayes and Empirical Bayes Methods 
# for Data Analysis, Chapman & Hall/CRC, pages 233-234.
#
################################################################################
g = function(gamma,x){x/(gamma+x)}
log.full.gamma = function(beta,gamma,sigma2){
  dnorm(gamma,g0,tau0,log=TRUE)+
  sum(dnorm(y,beta[1]+beta[2]*g(gamma,x),sqrt(sigma2),log=TRUE))
}

x = c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)
y = c(76,47,97,107,123,139,159,152,191,201,207,200)
n = length(x)

pdf(file="data.pdf",width=6,height=5)
plot(x,y,ylim=range(y,q1,q2),pch=16,xlab="Substrate concentration (ppm)",
     ylab="Velocity of enzymatic reaction (counts/min/min)",col=2)
dev.off()

# Prior hyperparameters
p     = 2
b0    = c(50,170)
B0    = diag(3,p)
g0    = 0.05
tau20 = 0.01
nu0   = 5
sigma20 = 126
tau0  = sqrt(tau20)

# Sufficient statistics and other constants
par1  = (nu0+n)/2
iB0   = solve(B0)
iB0b0 = iB0%*%b0
nu0sigma20 = nu0*sigma20


# MCMC setup and nitial values
burnin = 10000
M      = 10000
skip   = 10
niter  = burnin+skip*M
beta   = b0
gamma  = g0
sigma2 = sigma20

# Running algorithms 1 or 2
set.seed(12345)
sdg = 0.1 # standard deviation of RW-MH proposal
for (algorithm in 1:2){
draws  = matrix(0,niter,4)
for (iter in 1:niter){
  X   = cbind(1,g(gamma,x))
  
  # sampling from sigma form 
  # sigma|x,y,beta,gamma (full conditional or from 
  # sigma|x,y,gamma (marginal)
  B1 = solve(t(X)%*%X+iB0)
  b1 = B1%*%(t(X)%*%y+iB0b0)
  if (algorithm==1){
    par2 = (nu0sigma20+sum((y-X%*%beta)^2))/2
  }else{
    par2 = (nu0sigma20+sum((y-X%*%b1)*y)+t(b0-b1)%*%iB0b0)/2
  }
  sigma2 = 1/rgamma(1,par1,par2)

  # sampling beta|x,y,sigma2,gamma
  beta = b1 + sqrt(sigma2)*t(chol(B1))%*%rnorm(2)

  # sampling gamma|x,y,beta,sigma2
  gamma1 = rnorm(1,gamma,sdg)
  lognum = log.full.gamma(beta,gamma1,sigma2)
  logden = log.full.gamma(beta,gamma,sigma2)
  logaccept = min(0,lognum-logden)
  if (log(runif(1))<logaccept){
    gamma = gamma1  
  }
  
  # Storing the output
  draws[iter,] = c(beta,gamma,sqrt(sigma2))
}
draws = draws[seq(burnin+1,niter,by=skip),]

# Posterior predictives
N     = 50
xxx   = seq(min(x),max(x),length=N)
E     = matrix(0,M,N)
pred = matrix(0,M,N)
for (i in 1:M){
   E[i,]    = draws[i,1]+draws[i,2]*g(draws[i,3],xxx)
   pred[i,] = rnorm(N,E[i,],draws[i,4])
}
q  = t(apply(pred,2,quantile,c(0.05,0.5,0.95)))
qE = t(apply(E,2,quantile,c(0.05,0.5,0.95)))

if (algorithm==1){
  draws1 = draws
  q1 = q
  qE1 = qE
}else{
  draws1 = draws
  q2 = q
  qE2 = qE
}
	
}

#####################################################################################
# Comparing algorithms
#####################################################################################

# Comparing posterior approximations
names = c("beta0","beta1","gamma","sigma2")
pdf(file="posterior.pdf",width=12,height=9)
par(mfrow=c(3,4),mai=c(0.6,0.6,0.2,0.2))
for (i in 1:4){
  range = range(draws1[,i],draws2[,i])
  ts.plot(draws2[,i],xlab="Iteration",ylab="",main=names[i],ylim=range)
  lines(draws1[,i],col=2)
}
for (i in 1:4){
  acf1 = acf(draws1[,i],main="",plot=FALSE,lag.max=100)
  acf2 = acf(draws2[,i],main="",plot=FALSE,lag.max=100)
  plot(acf1$acf,type="l",col=2,xlab="Lag",ylab="ACF")
  lines(acf2$acf)
}
legend("topright",legend=c("[sigma2]+[beta|sigma2]","[sigma2|beta]+[beta|sigma2]"),col=1:2,lty=1,bty="n")
for (i in 1:4){
  f1 = density(draws1[,i])
  f2 = density(draws2[,i])
  xrange = range(draws1[,i],draws2[,i])	
  yrange = range(f1$y,f2$y)
  plot(f1$x,f1$y,col=2,xlab="",main="",xlim=xrange,ylim=yrange,type="l",ylab="Density")
  lines(f2$x,f2$y)
}
dev.off()


pdf(file="posteriorpredictive.pdf",width=9,height=9)
# Comparing posterior predictive approximations
par(mfrow=c(1,1),mai=rep(1,4))
plot(x,y,ylim=range(y,q1,q2),pch=16,xlab="Substrate concentration (in ppm)",
ylab="Velocity of an enzymatic reaction (in counts/min/min)",col=2)
lines(xxx,q1[,1],col=2,lty=2)
lines(xxx,q1[,3],col=2,lty=2)
lines(xxx,qE1[,2],col=2)
lines(xxx,qE1[,1],col=2,lty=3)
lines(xxx,qE1[,3],col=2,lty=3)
lines(xxx,q2[,1],lty=2)
lines(xxx,q2[,3],lty=2)
lines(xxx,qE2[,2])
lines(xxx,qE2[,1],lty=3)
lines(xxx,qE2[,3],lty=3)
legend("topleft",legend=c("[sigma2]+[beta|sigma2]","[sigma2|beta]+[beta|sigma2]"),col=1:2,lty=1,bty="n")
legend("bottomright",legend=c(
"E(ynew|xnew,x,y)","90% CI for E(ynew|xnew,x,y)","90% CI for p(ynew|xnew,x,y)"),
lty=1:3,bty="n")
dev.off()