## Setup M-H algorithm
p <- function(mu) {
  dt(mu, 5) * prod(dt(y, 5, mu))
}
mh <- function(y, N=10000, sigma=1) {
  mu <- numeric(N)
  for (i in 1:(N-1)) {
    proposal <- mu[i] + rnorm(1, sd=sigma)
    r <- p(proposal)/p(mu[i])
    accept <- rbinom(1, 1, min(1,r))
    mu[i+1] <- if (accept) proposal else mu[i]
  }
  cat("Acceptance rate:", mean(diff(mu)!=0), "\n")
  mu
}

## Slide 10
y <- c(-1, 1, 5)
xx <- seq(-5,5,len=100)
plot(xx, Vectorize(p)(xx), type="l", col="blue", lwd=3, xlab=~mu, ylab="Posterior density")
mu <- mh(y)
plot(mu[1:100], type="l", xlab="t", ylab=~mu, las=1)
plot(mu, type="l", xlab="t", ylab=~mu, las=1)

## Slide 11
y <- c(-1, 1, 5)+40
mu <- mh(y)
plot(mu[1:100], type="l", xlab="t", ylab=~mu, las=1)
plot(mu, type="l", xlab="t", ylab=~mu, las=1)

## Slide 15
mu <- mh(y, sigma=15)
plot(mu[1:100], type="l", xlab="t", ylab=~mu, las=1)
plot(mu, type="l", xlab="t", ylab=~mu, las=1)

## Slide 16
mu <- mh(y, sigma=200)
plot(mu[1:100], type="l", xlab="t", ylab=~mu, las=1)
plot(mu, type="l", xlab="t", ylab=~mu, las=1)

## Slide 222
res <- 99
x <- y <- seq(-3, 3, len=res)
Z <- matrix(NA, nrow=res, ncol=res)
Sigma <- matrix(c(1, 0.5, 0.5, 1), 2)
Omega <- solve(Sigma)
for (i in 1:res) {
  for (j in 1:res) {
    b <- c(x[i], y[j])
    Z[i,j] <- -crossprod(b, Omega)%*%b
  }
}
B <- matrix(NA, 5, 2)
B[1,] <- c(1,1)
for (i in 2:nrow(B)) {
  B[i,1] <- rnorm(1, Sigma[1,2]*B[i-1,2], sqrt(1 - Sigma[1,2]^2))
  B[i,2] <- rnorm(1, Sigma[1,2]*B[i,1], sqrt(1 - Sigma[1,2]^2))
}
contour(x, y, Z, drawlabels=FALSE, levels=c(-5, -1, -1e-1, -1e-2), lwd=2, col="gray", xlab=expression(beta[1]), ylab=expression(beta[2]))
points(B[1,1], B[1,2], pch=19, col="red")
for (i in 2:nrow(B)) {
  points(B[i,1], B[i-1,2], pch=19, col="red")
  lines(c(B[i-1,1], B[i,1]), c(B[i-1,2], B[i-1,2]), col="red", lwd=2)
  points(B[i,1], B[i,2], pch=19, col="red")
  lines(c(B[i,1], B[i,1]), c(B[i-1,2], B[i,2]), col="red", lwd=2)
}

## Slide 27
res <- 99
x <- y <- seq(-3, 3, len=res)
Z <- matrix(NA, nrow=res, ncol=res)
Sigma <- matrix(c(1, 0.98, 0.98, 1), 2)
Omega <- solve(Sigma)
for (i in 1:res) {
  for (j in 1:res) {
    b <- c(x[i], y[j])
    Z[i,j] <- -crossprod(b, Omega)%*%b
  }
}
B <- matrix(NA, 5, 2)
B[1,] <- c(1,1)
for (i in 2:nrow(B)) {
  B[i,1] <- rnorm(1, Sigma[1,2]*B[i-1,2], sqrt(1 - Sigma[1,2]^2))
  B[i,2] <- rnorm(1, Sigma[1,2]*B[i,1], sqrt(1 - Sigma[1,2]^2))
}
contour(x, y, Z, drawlabels=FALSE, levels=c(-5, -1, -1e-1, -1e-2), lwd=2, col="gray", xlab=expression(beta[1]), ylab=expression(beta[2]))
points(B[1,1], B[1,2], pch=19, col="red")
for (i in 2:nrow(B)) {
  points(B[i,1], B[i-1,2], pch=19, col="red")
  lines(c(B[i-1,1], B[i,1]), c(B[i-1,2], B[i-1,2]), col="red", lwd=2)
  points(B[i,1], B[i,2], pch=19, col="red")
  lines(c(B[i,1], B[i,1]), c(B[i-1,2], B[i,2]), col="red", lwd=2)
}