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

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

## mh.norm:
##
## function for independent oe RW-MH sampling from a bivariate
## normal distriobution with zero mean, unit variance
## and correlation parameter rho -- with starting
## values and RW variance parameter given.  when muq = NULL
## then RW is used, otherwise an independet MH is used

mh.norm <- function(T, rho=1, x0=0, y0=0, muq=NULL, s2q=1, plotit=FALSE)
  {

    ## plot the bivariate normal that we're sampling from
    if(plotit) {
      if(is.null(muq)) main <- "RW-MH from a bivariate normal"
      else main <- "Indep-MH from a bivariate normal"
      plot(ellipse(rho), xlim=c(-3,3), ylim=c(-3,3),
           type="l", col=2, lty=2, lwd=2, main=main,
           cex.main = 1.5, cex.lab = 1.5, cex.axis=1.5)
    }
    
    ## initialise to starting values
    x <- rep(x0, T)
    y <- rep(y0, T)

    if(plotit) points(x[1], y[1])
    
    for(t in 2:T) {

      ## symmetric proposal for x
      if(plotit) readline(paste("x sample ", t, ": ", sep=""))
      if(is.null(muq)) xp <- rnorm(1, x[t-1], sqrt(s2q))
      else xp <- rnorm(1, muq, sqrt(s2q))
      if(plotit) lines(c(x[t-1], xp), c(y[t-1], y[t-1]), col=3, lty=2)

  
      ## acceptance ratio for x proposal
      A <- dnorm(xp, rho*y[t-1], sqrt(1-rho^2))
      A <- A/dnorm(x[t-1], rho*y[t-1], sqrt(1-rho^2))
      alpha <- min(A, 1)

      ## accept or reject x proposal
      if(runif(1) < alpha) {
        x[t] <- xp
        if(plotit) lines(c(x[t-1], x[t]), c(y[t-1], y[t-1]))
      } else x[t] <- x[t-1]

      ## symmetric proposal for y
      if(plotit) readline(paste("y sample ", t, ": ", sep=""))
      if(is.null(muq)) yp <- rnorm(1, y[t-1], sqrt(s2q))
      else yp <- rnorm(1, muq, sqrt(s2q))
      if(plotit) lines(c(x[t], x[t]), c(y[t-1], yp), col=3, lty=2)

      ## acceptance ratio for y proposal
      A <- dnorm(yp, rho*x[t], sqrt(1-rho^2))
      A <- A/dnorm(y[t-1], rho*x[t], sqrt(1-rho^2))
      alpha <- min(A, 1)

      ## accept or reject y proposal
      if(runif(1) < alpha) {
        y[t] <- yp
        if(plotit) lines(c(x[t], x[t]), c(y[t-1], y[t]))
      } else  y[t] <- y[t-1]

       ## add the next sample to the plot
      if(plotit) points(x[t], y[t])
    }

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