## code to accompany the ACF & ESS example from
## lecture 10 of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy

## load the library with the effectiveSize function
## you may need to install.packages("coda")
library(coda)

## read in the Gibbs Sampling and MH code for the
## bivariate normal distribution
source("gibbs_norm.R")
source("mh_norm.R")

## use the effectiveSize funtion to look at the ESS
## for the marginal samples of x

rho <- 0.25
T <- 2000

## Gibbs sampling
xy.gibbs <- gibbs.norm(T, rho, 10, 10)

## RW-MH
xy.rwmh <- mh.norm(T, rho, 2, -2, NULL, s2q=0.25)

## Indep-MH
xy.indepmh <- mh.norm(T, rho, 2, -2, muq=0, s2q=2)

## plot ACF
par(mfrow=c(3,1), bty="n")
acf(xy.gibbs[,1], main="gibbs", lag.max=20)
acf(xy.rwmh[,1], main="RWM", lag.max=20)
acf(xy.indepmh[,1], main="Indep MH", lag.max=20)

## collect ESS into a data frame
ess.gibbs <- effectiveSize(xy.gibbs[,1])
ess.rwmh <- effectiveSize(xy.rwmh[,1])
ess.indepmh <- effectiveSize(xy.indepmh[,1])
ess <- data.frame(gibbs=ess.gibbs, rwmh=ess.rwmh,
                  indepmh=ess.indepmh)

## we can see that the ESS of the RW-MH algorithm
## and the Indep-MH algortithm is sensitive to
## the choice of sq2

## comparing RWMH with different choices of proposal s2q
T <- 1000
par(mfrow=c(3,2))
xy.001 <- mh.norm(T, rho, 2, -2, NULL, s2q=0.001)
plot(xy.001[,1], type="l", main="s2q=0.001", xlab="s", ylab="x", bty="n")
xy.01 <- mh.norm(T, rho, 2, -2, NULL, s2q=0.01)
plot(xy.01[,1], type="l", main="s2q=0.01", xlab="s", ylab="x", bty="n")
xy.1 <- mh.norm(T, rho, 2, -2, NULL, s2q=0.1)
plot(xy.1[,1], type="l", main="s2q=0.1", xlab="s", ylab="x", bty="n")
xy.10 <- mh.norm(T, rho, 2, -2, NULL, s2q=1)
plot(xy.10[,1], type="l", main="s2q=1", xlab="s", ylab="x", bty="n")
xy.100 <- mh.norm(T, rho, 2, -2, NULL, s2q=10)
plot(xy.100[,1], type="l", main="s2q=10", xlab="s", ylab="x", bty="n")
xy.1000 <- mh.norm(T, rho, 2, -2, NULL, s2q=100)
plot(xy.1000[,1], type="l", main="s2q=100", xlab="s", ylab="x", bty="n")