# Based on the examples in S Chib and E Greenberg, 
# Understanding the Metropolis-Hastings Algorithm,
# American Sytatistician 49 (1995), 327-335.
# 
library(lattice)
library(MASS)
n <- 600
Sigma <- matrix(c(1,0.9,0.9,1),nr=2)
InvSigma <- solve(Sigma)
P <- matrix(c(1,0,0.9,sqrt(1-0.9^2)),nr=2)
# 
# P can be used for the Choleski approach.
# Of course direct sampling using mvrnorm 
# (which presumably uses P) is possible.
# 
mu <- c(1,2)
lik <- function(x) as.numeric(exp(-(x-mu)%*%InvSigma%*%(x-mu)/2))
alpha <- function(x,y) min(lik(y)/lik(x),1)
xgrid <- -50:50/10
ygrid <- -50:50/10
mgrid <- matrix(0,length(xgrid),length(ygrid))
for (i in 1:101)
  for (j in 1:101){
    v <- cbind(xgrid[i],ygrid[j])
    mgrid[i,j] <- as.numeric(v %*% InvSigma %*% t(v))
  }
lev <- qchisq(0.99,2)
x1 <- matrix(0,n,2)
x1[1,] <- mu
k1 <- 0
for (i in 2:n){
  x <- x1[i-1,]
  z <- c(0.75,1) - c(1.5*runif(1),2*runif(1))
  y <- x + z
  k <- rbinom(1,1,alpha(x,y))
  k1 <- k1 + k
  x1[i,] <- x + k*(y-x)
}
cat("Method 1 has acceptance rate",100*k1/n,"%\n")
plot(x1,xlim=c(-4,6),ylim=c(-3,7))
par(new=T)
contour(x=xgrid+mu[1],y=ygrid+mu[2],mgrid,lev=lev,
        xlim=c(-4,6),ylim=c(-3,7),labels="99%")
x2 <- matrix(0,n,2)
x2[1,] <- mu
k2 <- k
for (i in 2:n){
  x <- x1[i-1,]
  z <- c(rnorm(1,sd=sqrt(0.6)),rnorm(1,sqrt(0.4)))
  y <- x + z
  k <- rbinom(1,1,alpha(x,y))
  x2[i,] <- x + k*(y-x)
  k2 <- k2 + k
}
cat("Method 2 has acceptance rate",100*k2/n,"%\n")
plot(x2,xlim=c(-4,6),ylim=c(-3,7))
par(new=T)
contour(x=xgrid+mu[1],y=ygrid+mu[2],mgrid,lev=lev,
        xlim=c(-4,6),ylim=c(-3,7),labels="99%")
x3 <- matrix(0,n,2)
x3[1,] <- mu
k3 <- 0
kr3 <- 0
D <- matrix(c(2,0,0,2),nr=2)
f <- function(x){
  as.numeric((2*pi)^(-1)*det(Sigma)^(-1/2)*exp((x-mu)%*%Sigma%*%(x-mu)))
}
h <- function(x){
  as.numeric((2*pi)^(-1)*det(D)^(-1/2)*exp((x-mu)%*%D%*%(x-mu)))
}
const <- 0.9
for (i in 2:n){
  cand <- FALSE
  while (!cand){
    kr3 <- kr3 + 1
    Z <- mvrnorm(1,mu,D)
    u <- runif(1)
    if (u <= f(Z)/(const*h(Z))){
      y <- Z
      cand <- TRUE
    }
  }
  kr3 <- kr3 - 1
  C1 <- as.numeric(f(x) < const*h(x))
  C2 <- as.numeric(f(y) < const*h(y))
  u <- runif(1)
  if (C1==1) Alpha <- 1
  if (C1==0 && C2==1) Alpha <- const*h(x)/f(x)
  if (C1==0 && C2==0) Alpha <- min(f(y)*h(x)/(f(x)*h(y)),1)
  k <- rbinom(1,1,Alpha)
  k3 <- k3 + k
  x3[i,] <- x + k*(y-x)
}
cat("Method 3 has acceptance rate",100*k3/n,"%\n")
cat("and uses an average of",kr3/n,"extra steps in the A-R stage\n")
plot(x3,xlim=c(-4,6),ylim=c(-3,7))
par(new=T)
contour(x=xgrid+mu[1],y=ygrid+mu[2],mgrid,lev=lev,
        xlim=c(-4,6),ylim=c(-3,7),labels="99%")
x4 <- matrix(0,n,2)
x4[1,] <- mu
k4 <- 0
for (i in 2:n){
  x <- x4[i-1,]
  z <- 1 - 2*c(runif(1),runif(1))
  y <- mu - (x - mu) + z
  k <- rbinom(1,1,alpha(x,y))
  k4 <- k4 + k
  x4[i,] <- x + k*(y-x)
}
cat("Method 4 has acceptance rate",100*k4/n,"%\n")
plot(x4,xlim=c(-4,6),ylim=c(-3,7))
par(new=T)
contour(x=xgrid+mu[1],y=ygrid+mu[2],mgrid,lev=lev,
        xlim=c(-4,6),ylim=c(-3,7),labels="99%")