library(coda)


## Random walk Metropolis Hastings
mcmc_rw <- function(M, sd){

  # store samples here
  xsamples <- rep(NA, M)

  # store acceptance rates for sets of iterations
  thin  <- 10
  acceptrates <- rep(NA, (M/thin))

  yes <- 0
  no <- 0

  # specify a starting value
  xsamples[1] <- -10

  # MH-Iteration
  for(k in 2:M){
    # value of the past iteration
    old <- xsamples[k-1]

    # propose a new value based on the old one
    proposal <- rnorm(1, mean=old, sd=sd)

    # compute acceptance ratio
    log.target.ratio <- dnorm(proposal, mean=0, sd=1, log=TRUE) - 
			dnorm(old, mean=0, sd=1, log=TRUE)

    # the proposal ratio is equal to 1 as we have 
    #a symmetric proposal distribution
    log.proposal.ratio <- 0

    # get the acceptance probability (on log scale)
    alpha <- log.target.ratio+log.proposal.ratio

    # accept-reject step
    if(log(runif(1)) <= alpha){	
      # accept the proposed value
      xsamples[k] <- proposal
      # increase counter of accepted values
      yes <- yes + 1
    }
    else{
      # stay with the old value
      xsamples[k] <- old
      no <- no + 1
    }

    # every 10 iterations, record the acceptance rate
    if(k%%thin == 0) 
      acceptrates[k/thin] <- yes/(yes+no)*100
  }

  # acceptance rate
  cat("The acceptance rate is: ", round(yes/(yes+no)*100,2), "%\n", sep="")
  return(list(xsamples=xsamples, acceptrates=acceptrates))
}

set.seed(03122010)
s1 <- mcmc_rw(10000, sd=0.24)
s2 <- mcmc_rw(10000, sd=2.4)
s3 <- mcmc_rw(10000, sd=24)

### Use coda to check convergence

mcmc1 = as.mcmc(s1[[1]])
mcmc2 = as.mcmc(s2[[1]])
mcmc3 = as.mcmc(s3[[1]])

plot(mcmc1)
x11()
plot(mcmc2)
x11()
plot(mcmc3)

mcmc1 = as.mcmc(s1[[1]][100:10000])
mcmc2 = as.mcmc(s2[[1]][100:10000])
mcmc3 = as.mcmc(s3[[1]][100:10000])

geweke.diag(mcmc1)
geweke.diag(mcmc2)
geweke.diag(mcmc3)

autocorr.plot(mcmc1)
autocorr.plot(mcmc2)
autocorr.plot(mcmc3)

source("ess.R")
ess(mcmc1, 9900)
ess(mcmc2, 9900)
ess(mcmc3, 9900)