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


## gibbs.norm:
##
## function for Gibbs sampling from a bivariate
## normal distriobution with zero mean, unit variance
## and correlation parameter rho -- with starting
## values given

gibbs.norm <- function(T, rho=1, x0=0, y0=0, plotit=FALSE)
  {
    
    ## plot the 95% contour of the normal we are
    ## sampling from
    if(plotit) {
      plot(ellipse(rho), type="l", col=2, lty=2, lwd=2,
           main="Gibbs sampling from a bivariate normal",
           cex.main = 1.5, cex.lab = 1.5, cex.axis=1.5)
    }
    
    ## initialise to starting values
    x <- rep(x0, T)
    y <- rep(y0, T)

    ## do the Gibbs sampling
    if(plotit) points(x[1], y[1])
    for(t in 2:T) {

      ## sample the x direction conditional on y
      if(plotit) readline(paste("x and y sample ", t, ": ", sep=""))
      x[t] <- rnorm(1, rho*y[t-1], sqrt(1-rho^2))
      if(plotit) lines(c(x[t-1], x[t]), c(y[t-1], y[t-1]))

      ## sample the y direction conditional on x 
      y[t] <- rnorm(1, rho*x[t], sqrt(1-rho^2))
      if(plotit) {
        lines(c(x[t], x[t]), c(y[t-1], y[t]))
        points(x[t], y[t])
      }
    }

    ## return a matrix
    return(cbind(x,y))
  }