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

## prior settings and data
s2 <- 1; tau2.0 <- 10; mu.0 <- 5
y <- c(9.37, 10.18, 9.16, 11.60, 10.33)
n <- length(y)
ybar <- mean(y)

## RWMH proposal variace
s2.q <- 2

## storage for samples
S <- 10000
mu <- rep(NA, S)
mu[1] <- 0

## RWMH loop
for(s in 2:S) {

  ## proposal
  mu.star <- rnorm(1, mu[s-1], sqrt(s2.q))

  ## (log) likelihood and prior of proposal
  llik.star <- sum(dnorm(y, mu.star, sqrt(s2), log=TRUE))
  lprior.star <- dnorm(mu.star, mu.0, sqrt(tau2.0), log=TRUE)

  ## (log) likelihood and prior of previous value
  llik <- sum(dnorm(y, mu[s-1], sqrt(s2), log=TRUE))
  lprior <- dnorm(mu[s-1], mu.0, sqrt(tau2.0), log=TRUE)

  ## calculate the acceptance ratio
  A <- exp(llik.star + lprior.star - llik - lprior)

  ## accept or reject
  if(A > 1 || runif(1) < A) mu[s] <- mu.star ## accept
  else mu[s] <- mu[s-1] ## reject
}


## plot the chain
plot(mu, type="l", bty="n", ylab="mu", xlab="s")

## plot a histogram of the samples and compare to the truth
hist(mu[-(1:100)], bty="n", main="", xlab="mu", prob=TRUE)
o <- order(mu)

## calculate the true distribution and comare
tau2.n <- 1/(n/s2 + 1/tau2.0)
mu.n <- ybar*(n/s2)*tau2.n + mu.0*(1/tau2.0)*tau2.n
lines(mu[o], dnorm(mu[o], mu.n, sqrt(tau2.n)),
      col=2, lwd=2, lty=2)