Partial example: The Metropolis-Hastings sampler

Below, you can find an example of how to code a Metropolis-Hastings sampler. Some bits are left out for you to fill in (marked “INSERT HERE”). Each “INSERT HERE” statement requires one line of code. If you struggle, you can find a link to a solution below the function.

# This is a function that takes four parameters:
# - target: the target distribution, a function that takes one
#   argument (a number) and returns the (logged) value of a
#   distribution
# - initTheta: the initial value of theta, a number
# - proposalSd: the standard deviation of (Gaussian) proposal
#   distribution
# - nIterations: the number of iterations
# The function returns a vector of samples of theta from the target
# distribution
my_mcmcMh <- function(target, initTheta, proposalSd, nIterations) {
  # evaluate the function "target" at "initTheta", and assign to
  # a variable called targetThetaCurrent.
  targetThetaCurrent <- target(initTheta)

  # initialise variables to store the current value of theta, the
  # vector of samples, and the number of accepted proposals
  thetaCurrent <- initTheta
  samples <- initTheta
  accepted <- 0

  # run MCMC for nIteration interations
  for (i in seq_len(nIterations)) {
    thetaProposed <- # INSERT HERE: draw a new theta from the
      # (Gaussian) proposal distribution and assign to a
      # variable called "thetaProposed". See "?rnorm for help.

      targetThetaProposed <- # INSERT HERE: evaluate the function
      # target at the proposed theta and assign to a
      # variable called "targetThetaProposed"

      logAcceptance <- # INSERT HERE: compute Metropolis-Hastings ratio
      # (acceptance probability). This is easiest if you assume
      # the target function to return the logarithm of the
      # distribution value. Assign the result to a variable
      # called "logAcceptance"

      r <- # INSERT HERE: draw random number number between 0 and 1
      # using "runif" and assign to a variable called "r".

      # test acceptance by comparing the random number to the
      # Metropolis-Hastings ratio (acceptance probability) (using
      # "exp" because we calculated the logarithm of the
      # Metropolis-Hastings ratio before)
      if (r < exp(logAcceptance)) {
        # if accepted:
        thetaCurrent <- # INSERT HERE: change the current value
          # of theta to the proposed theta

          targetThetaCurrent <- # INSERT HERE: updated the current
          # value of the target

          # update number of accepted proposals
          accepted <- accepted + 1
      }

    # add the current theta to the vector of samples
    samples <- c(samples, thetaCurrent)

    # print current state of chain and acceptance rate
    cat(
      "chain:", thetaCurrent,
      "acceptance rate:", accepted / i, "\n"
    )
  }

  return(samples)
}

If you run into any problems, have a look at our solution. Otherwise you can return to the practical.

 

---

This web site and the material contained in it were originally created in support of an annual short course on Model Fitting and Inference for Infectious Disease Dynamics at the London School of Hygiene & Tropical Medicine. All material is under a MIT license. Please report any issues or suggestions for improvement on the corresponding GitHub issue tracker. We are always keen to hear about any uses of the material here, so please do get in touch using the Discussion board if you have any questions or ideas, or if you find the material here useful or use it in your own teaching.