#
# MCMC Diagnostics examples
# For STAT8810
# Fall, 2017
# M.T. Pratola



# Generate our usual GP example
source("dace.sim.r")
set.seed(88)
n=10
k=1
rhotrue=0.2
lambdaytrue=1
design=as.matrix(runif(n))
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)




# Run the MCMC to fit the Bayesian GP regression 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)

# regress works as follows:
# Adapt for first 50% of iterations -- here that is 2500.
# Further burn-in is draws up to start of last 
#      iterations -- here that is 2501--4000.
# last is number of draws to take as posterior samples.  
#      Here that is 4001--5000.

fit=regress(y,l.dez,5000,pi,mh,last=1000,adapt=FALSE)





# Trace plots
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 the MCMC with a proposal width of 1e-6, no adaptation.
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)

# Run the MCMC
fit=regress(y,l.dez,20000,pi,mh,last=5000,adapt=FALSE)




# Traceplots
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)




# Autocorrelation plots
par(mfrow=c(1,2))
acf(fit$lambdaz,main="lambdaz")
acf(fit$rhoz,main="rhoz")





# Run the MCMC again, but here I will return everything.
fit=regress(y,l.dez,20000,pi,mh,last=20000,adapt=FALSE)


# Raftery & Lewis
library(coda)
post=as.matrix(cbind(fit$lambdaz,fit$rhoz))
postmcmc=as.mcmc(post)
raftery.diag(postmcmc,q=0.025)
raftery.diag(postmcmc,q=0.975)





# Run the MCMC, but now turn on adaptation.  Again,
# we will return everything.
fit2=regress(y,l.dez,20000,pi,mh,last=20000,adapt=TRUE)

# Raftery & Lewis
library(coda)
post=as.matrix(cbind(fit2$lambdaz,fit2$rhoz))
postmcmc=as.mcmc(post)

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


# ACF plots
par(mfrow=c(1,2))
acf(fit2$lambdaz,main="lambdaz")
acf(fit2$rhoz,main="rhoz")


# R&L diagnostic
raftery.diag(postmcmc,q=0.025)
raftery.diag(postmcmc,q=0.975)





# Run the MCMC, but now turn on adaptation. 
# Now just return the last 5000 since it looks like
# we burn-in by then.
fit3=regress(y,l.dez,20000,pi,mh,last=5000,adapt=TRUE)

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

# ACF plots
par(mfrow=c(1,2))
acf(fit3$lambdaz,main="lambdaz")
acf(fit3$rhoz,main="rhoz")

# Raftery & Lewis
library(coda)
post=as.matrix(cbind(fit3$lambdaz,fit3$rhoz))
postmcmc=as.mcmc(post)
raftery.diag(postmcmc,q=0.025)
raftery.diag(postmcmc,q=0.975)




# Geweke Diagnostic
post=as.matrix(cbind(fit$lambdaz,fit$rhoz))
postmcmc=as.mcmc(post)
geweke.diag(post)

post=as.matrix(cbind(fit2$lambdaz,fit2$rhoz))
postmcmc=as.mcmc(post)
geweke.diag(post)

post=as.matrix(cbind(fit3$lambdaz,fit3$rhoz))
postmcmc=as.mcmc(post)
geweke.diag(post)





# Gelman-Rubin Diagnostic
postdraws=vector("list",6)
fit=vector("list",6)
for(i in 1:6) fit[[i]]=regress(y,l.dez,20000,pi,mh,last=5000,adapt=TRUE)
for(i in 1:6) postdraws[[i]]=fit[[i]]$rhoz
for(i in 1:6) postdraws[[i]]=as.mcmc(postrows[[i]])
postmulti=as.mcmc.list(postrows)
gelman.diag(postmulti,autoburnin=FALSE)





# Effective Sample Size
post=as.matrix(cbind(fit2$lambdaz,fit2$rhoz))
postmcmc=as.mcmc(post)
effectiveSize(postmcmc)

post=as.matrix(cbind(fit3$lambdaz,fit3$rhoz))
postmcmc=as.mcmc(post)
effectiveSize(postmcmc)