library(coda) #Install package coda!

############################################################
## Markov Chain Monte Carlo
## CO902 University of Warwick 2014 - Ben Graham
############################################################
## Code to run a Markov chain 
## x, iterative.function(x), iterative.function(iterative.function(x)), ...
## Allows a burn-in period to reduce bias caused by starting 
## point. Allows thinning to reduce memory when slow mixing 
## forces you to choose a large number of iterations
############################################################
run.mcmc <- function(
  #Stochastic map from x[i] to x[i+1]
  iterative.function,
  #start point
  x, 
  #how long to run for
  iterations=1000,
  #discard initial output to allow Markov chain mixing to take place
  burn.in=0, 
  #thinning factor, >=1, store a subset of the output of slow mixing chains to save memory
  thin=1) {  
  ##run iterative.function for n steps, starting from x
  steps <- function(x,n) {
    if (n>0)
      for (i in 1:n)
        x=iterative.function(x)
    x
  }
  n=floor((iterations-burn.in)/thin)
  #Place to store output
  chain=matrix(0,n,length(x))
  #Burn-in the chain
  x=steps(x,burn.in)
  chain[1,]=x
  if (n>2)
    for (i in 2:n)
      chain[i,]=steps(chain[i-1,],thin)
  #Return an coda mcmc object
  mcmc(chain,start=burn.in,thin=thin)
}

metropolis.hastings.mcmc.sampler <- function(
  #log target density
  log.pi,
  #random function: given x, propose a y
  rproposal,
  #log density of the proposal distribution
  log.dproposal,
  ...) {
  iterative.function <- function(x) {
    y=rproposal(x)
    alpha=min(1, exp(log.pi(y)+log.dproposal(y,x)
                     -log.pi(x)-log.dproposal(x,y)))
    if (runif(1)<alpha) y
    else x
  }
  function(start.point) {
  run.mcmc(iterative.function,
           start.point,...)
  }
}

independence.sampler.markov.chain <- function(
  #target density
  log.pi,
  #random function that does not take any arguments
  rproposal,
  #log density of the proposal distribution
  log.dproposal,...) {
    metropolis.hastings.mcmc.sampler(
    log.pi,
    function(x) rproposal(),
    function(x,y) log.dproposal(y),...)
}


random.walk.sampler.markov.chain <- function(
  #target density
  log.pi,
  #Proposal distribution; Take N(x,sd^2) steps 
  sd,...) {
  metropolis.hastings.mcmc.sampler(
    log.pi,
    function(x) x+rnorm(length(x),sd=sd),
    function(x,y) sum(log(dnorm(y-x,sd=sd))),...)
}

gibbs.sampler.markov.chain <- function(
  target.conditional.distributions,...) {
  iterative.function <- function(x) {
    i=sample(1:length(x),1)
    x[i]=(target.conditional.distributions[[i]])(x)
    x
  }
  function(start.point)
    run.mcmc(iterative.function,
             start.point,...)
}

run.gelman.diagnostic <- function(x.list,fn) {
  l=mcmc.list(lapply(x.list,fn))
  print(gelman.diag(l,multivariate = FALSE))
  gelman.plot(l)
}

