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

## load the ellipse package
## you will need to install.packages("ellipse")
library(ellipse)

## read the file with the gibbs.norm function in it
source("gibbs_norm.R")

## want T=50 samples from the bivariate
## normal distribution with unit variance,
## zero mean, and correlation rho
rho <- 0.25
xy <- gibbs.norm(50, rho, plotit=TRUE)

## then takeT=2000 samples starting at (10,10)
T <- 2000
xy <- gibbs.norm(T, rho, 10, 10)

## plot the samples
plot(xy, main="Gibbs sampling 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[10:T,2])

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