## code to accompany the MH sampling example from
## lecture 9 of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy


## Using Metropolis Hastings (MH) to sample
## from a bivariate normal (posterior) distribution
source("mh_norm.R")

## want to sample from the bivariate
## normal distribution with unit variance,
## zero mean, and correlation rho
rho <- 0.25

##
## start with the Randowm Walk (RW)-MH using a
## symmetric normal proposal centered at the
## previous value with variance s2q
##

s2q <- 0.25

## illustrate by stepping through a few rounds
xy <- mh.norm(50, rho, 2, -2, NULL, s2q, TRUE)

## show the traces of the dependent samples
par(mfrow=c(2,1), cex.main=1.5, cex.lab=1.5,
    cex.axis=1.5)
plot(xy[,1], type="b", main="RW-MH samples of x")
plot(xy[,2], type="b", main="RW-MH samples of y")

## then take T=2000 samples starting at (2,-2)
T <- 2000
xy <- mh.norm(T, rho, 2, -2, NULL, s2q)

## plot the samples
par(mfrow=c(1,1))
plot(xy[100:T,], main="RW-MH from a bivariate normal",
     cex.main = 1.5, cex.lab = 1.5, cex.axis=1.5)
lines(ellipse(rho), col=2, lty=2, lwd=2)

## look at a plot of the trace of the x and y
## and put histograms on the right hand side
par(mfrow=c(2,2), cex.main=1.5, cex.axis=1.5,
    cex.lab=1.5)

## trace and histogram for x
plot(xy[,1], type="l", main="trace of x")
hist(xy[100:T,1])

## trace and histogram for y
plot(xy[,2], type="l", main="trace of y")
hist(xy[100:T,2])

## calculate the correlation parameter
cor(xy[100:T,])

##
## Now use an independence MH with normal proposal
## centered at zero with variance 2
##

par(mfrow=c(1,1), cex.main=1.5, cex.axis=1.5,
    cex.lab=1.5)
muq <- 0
s2q <- 2

## illustrate by stepping through a few rounds
xy <- mh.norm(50, rho, 2, -2, muq, s2q, TRUE)

## show the traces of the dependent samples
par(mfrow=c(2,1), cex.main=1.5, cex.lab=1.5,
    cex.axis=1.5)
plot(xy[,1], type="b", main="Indep-MH samples of x")
plot(xy[,2], type="b", main="Indep-MH samples of y")

## then take T=2000 samples starting at (2,-2)
T <- 2000
xy <- mh.norm(T, rho, 2, -2, muq, s2q)

## plot the samples
par(mfrow=c(1,1))
plot(xy, main="Indep-MH from a bivariate normal",
     cex.main = 1.5, cex.lab = 1.5, cex.axis=1.5)
lines(ellipse(rho), col=2, lty=2, lwd=2)

## look at a plot of the trace of the x and y
## and put histograms on the right hand side
par(mfrow=c(2,2), cex.main=1.5, cex.axis=1.5,
    cex.lab=1.5)

## trace and histogram for x
plot(xy[,1], type="l", main="trace of x")
hist(xy[10:T,1])

## trace and histogram for y
plot(xy[,2], type="l", main="trace of y")
hist(xy[100:T,2])

## calculate the correlation parameter
cor(xy[100:T,])


## LATER, return to this example and
## talk about effective sample size