#
# Bayesian Gaussian Process Example
# For STAT8810
# Fall, 2017
# M.T. Pratola


##########################################################################################################################
# First, generate the response.
##########################################################################################################################
cat("\nGenerating response\n")
n=10
k=1
rhotrue=0.2
lambdaytrue=1
design=as.matrix(runif.sobol(n,1))
l1=list(m1=outer(design[,1],design[,1],"-"))
l.dez=list(l1=l1)
R=rhogeodacecormat(l.dez,c(rhotrue))$R
L=chol(R)
u=rnorm(nrow(R))
z=t(L)%*%u

# now set up our observed data:
y=z
l1=list(m1=outer(design[,1],design[,1],"-"))
l.dez=list(l1=l1)


##########################################################################################################################
# Now fit the GP model
##########################################################################################################################
source("regression.r")
pi=list(az=5,bz=5,rhoa=rep(1,k),rhob=rep(5,k))


# Run MCMC using a proposal width of 1e-5 for the rho parameters.
mh=list(rr=1e-5)
fit=regress(y,l.dez,5000,pi,mh,last=2000,adapt=FALSE)

par(mfrow=c(1,2))
plot(fit$lambdaz,type='l',xlab="draw",ylab="lambdaz")
abline(h=lambdaytrue)
plot(fit$rhoz,type='l',xlab="draw",ylab="rhoz")
abline(h=rhotrue)


# Run MCMC using a proposal width of 0.5 for the rho parameters.
mh=list(rr=.5)
fit=regress(y,l.dez,5000,pi,mh,last=2000,adapt=FALSE)

par(mfrow=c(1,2))
plot(fit$lambdaz,type='l',xlab="draw",ylab="lambdaz",ylim=c(0,2.5))
abline(h=lambdaytrue)
plot(fit$rhoz,type='l',xlab="draw",ylab="rhoz",ylim=c(0,1))
abline(h=rhotrue)

# Run MCMC using a proposal width of 0.27 for the rho parameters.
mh=list(rr=0.27)
fit=regress(y,l.dez,5000,pi,mh,last=2000,adapt=FALSE)

par(mfrow=c(1,2))
plot(fit$lambdaz,type='l',xlab="draw",ylab="lambdaz",ylim=c(0,2.5))
abline(h=lambdaytrue)
plot(fit$rhoz,type='l',xlab="draw",ylab="rhoz",ylim=c(0,1))
abline(h=rhotrue)

# Run MCMC using a adaptive MCMC with a starting proposal width of 0.05 for the rho parameters.
mh=list(rr=.05)
fit=regress(y,l.dez,5000,pi,mh,last=2000,adapt=TRUE)

par(mfrow=c(1,2))
plot(fit$lambdaz,type='l',xlab="draw",ylab="lambdaz",ylim=c(0,2.5))
abline(h=lambdaytrue)
plot(fit$rhoz,type='l',xlab="draw",ylab="rhoz",ylim=c(0,1))
abline(h=rhotrue)




# Predict the last fit on a grid of 100 equally spaced points
grid=as.matrix(seq(0,1,length=100))
design.all=rbind(design,grid)
l1=list(m1=outer(design.all[,1],design.all[,1],"-"))
l.v=list(l1=l1)
fitp=predict(l.v,fit)

# Plot the posterior predictions using mean and sd to quantify uncertainty (pointwise)
plot(design,y,pch=20,col="red",cex=2,xlim=c(0,1),
    ylim=c(2.5,3.5),xlab="x",
    main="Predicted mean response +/- 2s.d.")
for(i in 1:nrow(fitp$preds))
  lines(grid,fitp$preds[i,],col="grey",lwd=0.25)
mean=apply(fitp$preds,2,mean)
sd=apply(fitp$preds,2,sd)
lines(grid,mean-1.96*sd,lwd=0.75,col="black")
lines(grid,mean+1.96*sd,lwd=0.75,col="black")
lines(grid,mean,lwd=2,col="blue")

# Plot the posterior predictions using median and quantiles to quantify uncertainty (pointwise)
plot(design,y,pch=20,col="red",cex=2,xlim=c(0,1),ylim=c(2.5,3.5),
    xlab="x",main="Predicted median, q.025 and q.975")
for(i in 1:nrow(fitp$preds))
  lines(grid,fitp$preds[i,],col="grey",lwd=0.25)
med=apply(fitp$preds,2,quantile,0.5)
q.025=apply(fitp$preds,2,quantile,0.025)
q.975=apply(fitp$preds,2,quantile,0.975)
lines(grid,q.025,lwd=0.75,col="black")
lines(grid,q.975,lwd=0.75,col="black")
lines(grid,mean,lwd=2,col="blue")