require(coda)

## a fairly general Metropolis-Hastings sampler
## this code assumes that the proposal distribution is symmetric

mcmc.mh <- function (theta, fun, niter, thin = 1, burnin = 0, progress = FALSE, rprop = NULL) {
  if (inherits(theta,'mcmc')) {
    nm <- colnames(theta)
    en <- end(theta)
    theta <- theta[nrow(theta),]
    names(theta) <- nm
  } else {
    en <- 0
  }
  if (is.null(rprop))
    rprop <- function(x)runif(n=length(x),min=x-10,max=x+10)
  chain <- matrix(
                  nrow=as.integer(niter),
                  ncol=length(theta),
                  dimnames=list(NULL,names(theta))
                  )
  f <- fun(theta)
  k <- 1                                # step in the chain
  nsave <- 0                            # number of values saved so far
  nprop <- 0                            # number of proposals made
  nacpt <- 0                            # number of proposals accepted
  if (progress) {                       # do some progress-bar-related stuff
    pb <- txtProgressBar(title="mcmc.mh",label="progress",style=3)
    npb <- floor((burnin+niter)/100)
  }
  while (nsave < niter) {
    theta.star <- rprop(theta)          # propose a new parameter
    f.star <- fun(theta.star)
    nprop <- nprop+1
    alpha <- f.star/f  # Metropolis' rule (since proposal is symmetric)
    if ((alpha>=1)||(runif(1)<alpha)) {
      theta <- theta.star
      f <- f.star
      nacpt <- nacpt+1
    }
    k <- k+1
    if ((k>burnin)&&(k%%thin==0)) {
      nsave <- nsave+1
      chain[nsave,] <- theta
    }
    if (progress&&(k%%npb==0))
      setTxtProgressBar(pb,value=k/(niter*thin+burnin))
  }
  if (progress)
    close(pb)
  cat("acceptance fraction = ",nacpt/nprop,"\n")
  mcmc(data=chain,thin=thin,start=burnin+en+1) # store the result in a CODA object
}