Why in Hamiltonian MCMC do we multiply the posterior distribution by the likelihood?

Question

So maybe I am misunderstanding what the author is staying, but I am reading Chapter 14 of Kruschke's Doing Bayesian Analysis. I am reading about the software Stan and how it uses the Hamiltonian MCMC (HMC) to sample from the posterior distribution. To recap, HMC generates a proposal (jumping) distribution based on the gradient of the posterior distribution. To do HMC in Stan, we specify a model by defining both a distribution (or likelihood) for the data and the prior, for example:

model {
    theta ~ beta(1,1) ;  // prior
    y ~ bernoulli(theta) // likelihood
  }

But then I read a blurb that says:

In other words, in Stan, we are not directly randomly sampling from the posterior (like we would do in Gibbs where we sample from the conditional posterior). Instead, we determine a proposal point via the gradient of the posterior. What I don't get is why are we multiplying the posterior by the likelihood? I don't understand.

Thanks!

Answer 1

Stan computes a log-posterior density and uses its gradient to do sampling. It does this by incrementing a variable storing the log probability (really, the log kernel. Ben Goodrich points out that Stan only needs to care about the log probability up to constant terms, which are neglected). At each iteration, each sampling statement in the model block increments the log probability variable, so when you have several sampling statements, the value of the log posterior is incremented several times. When the author writes "current posterior" in the quote, the author means "the current value of the log posterior accumulator." That's what the "running total" in the last sentence refers to.

So to answer your question

[W]hy are we multiplying the posterior by the likelihood?

Instead of a "running total" on the log scale, we can do equivalent computations using a "running product" on the original scale, where we compute the posterior density in stages, one product at a time.

Stan still does the computations "under the hood" by incrementing a log probability, but as a matter of convenience, users can write y ~ distribution(arguments) instead of writing increment_log_prob(distribution_log(arguments)) .