The Proposal Step

The proposal step starts with a initial guess at \(\mu\), which we call \(\mu^{(0)}\). We are going to have some fun and initialize poorly. Let’s set \(\mu^{(0)} = .1\).

Once we have our initial sample, we need to update that guess with a proposal. We propose a new value \(\mu^{*}\) from a proposal distribution:

\[\mu^{*}|\mu^{(0)} \sim N( \mu^{(0)}, \tau^2)\]

For now, let’s set \(\tau^2 = (.2^2) =.04\). We’ll explore other choices of \(\tau\) later.

When a proposal is not valid, we discard the proposal and try again.

The Accept/Reject Step

Now that we have a proposal of \(\mu^{*}\), we have a decision to make. Do we think that \(\mu^{*}\) is more likely to be draw from the posterior of \(\mu\) than \(\mu^{(0)}\)? If so, we move our chain to \(\mu^{*}\) by setting

\[\mu^{(1)} = \mu^{*}\]

However, if \(\mu^{*}\) is not likely draw from the posterior, than we want to leave our chain right where it is by setting

\[\mu^{(1)} = \mu^{(0)}\]

This means we somehow need to determine whether or not \(\mu^{*}\) is likely to be a draw from the posterior (or at least, more likely than \(\mu^{(0)}\)).

Typically, we don’t have the actual posterior (if we did, we wouldn’t need Metropolis). In Metropolis, we use an acceptance ratio to help us determine the next step in our MCMC chain.

\[R = min\left( \frac{P(\mu^{*} | Y_1, \dots, Y_n)}{P(\mu^{(0)} | Y_1, \dots, Y_n)}, 1 \right) = min( R_{inner}, 1)\] We can simplify the inner part of the ratio, \(R_{inner}\), to

\[\begin{aligned} R_{inner} &= \frac{P(Y_1, \dots, Y_n|\mu^{*})P(\mu^{*} )}{P(Y_1, \dots, Y_n|\mu^{(0)})P(\mu^{(0)} )} \\ &= \underbrace{\frac{P(Y_1, \dots, Y_n|\mu^{*})}{P(Y_1, \dots, Y_n|\mu^{(0)})}}_\text{Likelihood Component} \times \overbrace{\frac{P(\mu^{*} )}{P(\mu^{(0)} )}}^\text{Prior Component} \end{aligned}\]

With both components computed, we now know the form of the acceptance ratio!! The next step is to write a function that will help us to compute it.

Pause: What’s a function?

A function is a computer structure that takes inputs and produces an output. It has a structure that looks like this:

nameOfFunction <- function( inputs ){
  
  Compute what you want to output 
  
  return(output)

}

In our case, we want to input \(\mu^{*}\) and \(\mu = \mu^{(0)}\) and output our acceptance ratio R. We also need to input any constants, like \(a\) or \(\sigma\), that are needed in the ratio. This means we will have a structure like

computeAcceptanceRatio <- function( mustar, mu , a, b, sigma, n, ybar ){
  
  Compute R 
  
  return(R)
}

Inside the function we need to put in the steps to compute R, which luckily we just derived!

computeAcceptanceRatio <- function( mustar, mu , a, b, sigma, n, ybar ){
  
  # Compute the likelihood part
  likelihood_numerator <- 
  likelihood_denominator <- 
  likelihood_component <- likelihood_numerator / likelihood_denominator
  
  # Compute the prior part
  prior_numerator <- 
  prior_denominator <- 
  prior_component <- prior_numerator / prior_denominator
  
  # Multiply the two
  R <- likelihood_component * prior_component
  
  # Store 1 or the value from previous step
  R <- min( 1, R)
  
  # Output R
  return(R)
}

The value \(R\) you have just computed is the probability of accepting \(\mu^{*}\). This means that \((R \times 100)\)% of the time we propose \(\mu^{*}\) we want to accept it, and \(((1-R)\times 100)\)% of the time we want to reject it. How can we do that??

We are going to draw a random number \(u\) from a \(Uniform(0,1)\). If \(u\) is less than \(R\), then we accept \(\mu^{*}\). Otherwise, we reject \(\mu^{*}\)

And that’s one iteration of Metropolis!

Turning in your assignment

When your Markdown document is complete, do a final knit and make sure everything compiles. You will then submit your document on Canvas. You must submit a PDF or html document to Canvas. No other formats will be accepted. Make sure you look through the final document to make sure everything has knit correctly.