12.3 Concluding points

While we have done Metropolis-Hastings “by hand,” you may realize that we can automate this process. We will explore that in the next section. However there are several modifications we can do to make this algorithm more sophisticated:

• While we have focused on implementation of the Metropolis-Hastings algorithm with one parameter, this is easily extended to sets of parameter values, but we just change one parameter at a time.
• We can select parameter values from multiple starting points to make sure we don’t happen to start a chain near a local optimum. These starting points are called chains, and we can select the chain that has the highest likelihood value.
• We can use log likelihoods to make the likelihood ratio computation easier. By this way, instead of division we compute a difference of log likelihoods. This becomes numerically easier for R to handle, especially when the likelihood values are near zero. Here is some sample code that you can do this:

# Define a function that computes the LOG likelihood ratio for Wilson's weight
log_likelihood_ratio_wilson <- function(proposed, current) {

  # Define the model we are using
  my_model <- mass ~ p1 / (1 + exp(p2 - p3 * days))


  # This allows for all the possible combinations of parameters
  parameters <- tibble(p1 = c(current, proposed), p2 = 2.461935, p3 = 0.017032)

  # Compute the likelihood and return the likelihood from the list
  # Notice we've set logLikely = TRUE to compute the log likelihood
  out_values <- compute_likelihood(my_model, wilson, parameters, logLikely = TRUE)$likelihood

  # Return the likelihood from the list - here we compute the DIFFERENCE of likelihoods:
  ll_ratio <- out_values$l_hood[[2]] - out_values$l_hood[[1]]

  return(ll_ratio)
}

• As you would expect, the more times we iterate through this process, the better. Your initial guesses probably weren’t that great (or close to the global optimum), so a common procedure is to throw out the first percentage of iterations and call that the “burn-in” period.
• We can systematically explore the parameter space, where the jump distance changes depending on if we are always accepting new parameters or not. This process has several different implementations, but one is called simulated annealing.