I'm trying to understand how to estimate the parameter vector $\mathbf{\theta} = (\theta_1,\theta_2, \theta_3)$ of a model using the MH algorithm.

I am given a joint posterior density:

$p(\mathbf{\theta}, \alpha | y) \propto \alpha^n e^{\left(-\alpha f(\mathbf{\theta})\right)}$

I am also given two conditional posterior distributions:

  1. $p_{1}(\alpha|\theta,y) \propto IG(n,f(\theta))$
  2. $p_{2}(\theta|y) \propto f(\theta)^{-n}$

Where $f(\theta)$ is always the same function (cannot sample from it directly) which depends on $(\theta_1, \theta_2, \theta_3)$. Also, the parameter $\alpha$ is of no interest.

This is how I'm implementing the algorithm (the proposal $q$ is a multivariate Gaussian distribution):

$\theta^* \sim \mathcal{N}(\theta^{i-1}, \Sigma)$

$p = \min \left( \frac{p_2(\theta^*|y) q(\theta^{i-1}|\theta^*)}{p_2(\theta^{i-1}|y) q(\theta^*|\theta^{i-1})},1 \right) $

accept $\theta^*$ with probability $p$.

Where $q(\theta^{i-1}|\theta^*)$ is the density of the multivariate Gaussian distribution centered at $\theta^*$.

Question : am I supposed to be using the conditional posterior density $p_2$ for this step or should I use the joint posterior density? What is the logic behind this? So far the book examples I've seen use the joint posterior, but this is my first try with conditional posteriors (in non standard form) so I'm not sure I understand the logic.

Question 2 : should I be sampling from my proposal all at once or should I make a decision for each individual $\theta_i$?

Question 3 : should I be using the joint posterior somewhere ? Or is it needed to obtain the conditional posteriors only?

  • $\begingroup$ Q1: In our book we suggest simulation from both conditional posteriors because the marginal posterior is not often available. When it is, simulating directly from $p_2$ is more efficient. $\endgroup$
    – Xi'an
    Jun 22 '18 at 20:38
  • $\begingroup$ Q2: one can use a Gibbs sampler on the components of $\theta$ or a joint proposal, depending how easy it is. $\endgroup$
    – Xi'an
    Jun 22 '18 at 20:40

Question 1: Note that $p_2(\theta|y)$ is obtained by integrating out $\alpha$ from $p(\theta, \alpha|y)$. Since you don't care about $\alpha$, this is what you should use - it's the marginal distribution of the vector $\theta | y$ (not the conditional distribution, if by "conditional" you are referring to $\alpha$ rather than $y$.) If you sample from the joint posterior $p(\theta, \alpha|y)$, then ignore $\alpha$, you'll be integrating out $\alpha$ by Monte Carlo instead of by calculus, and, as one might expect, this won't be as accurate.

Question 2: You should sample from your proposal all at once UNLESS you are making separate proposals for each element of $\theta$. If you are making one proposal for $\theta$, the calculations of $q$ will include the effects of the proposals for all the elements of $\theta$, which will break the MCMC sampler's convergence to the posterior if you are actually only going to update one element of $\theta$.

The choice between making proposals for each element of $\theta$ vs. making one joint proposal is really an empirical one that is problem-specific. A block update that has a reasonable average acceptance probability (e.g., 15-40% across the entire run) can be more efficient at covering the region of the parameter space where the posterior lies than separate updates of each element, especially if the elements are correlated and the proposal distribution reflects some or most of that correlation; however, it is often harder to tune the proposal distribution to achieve reasonable acceptance probabilities if it is multivariate rather than a collection of univariate ones.

Question 3: See question 1; it was only needed to obtain the marginal distribution of $\theta|y$ (marginal with respect to $\alpha$, that is.)

  • $\begingroup$ Actually simulating from $p(\alpha,\theta|y)$ or from $p(\theta|y)$ produces exactly the same outcome in $\theta$: generations from $p(\theta|y)$. $\endgroup$
    – Xi'an
    Jul 21 '18 at 7:36
  • $\begingroup$ @Xi'an - yes, you're just integrating out $\alpha$ by Monte Carlo, but insofar as that integration isn't perfect for any finite sample size, your estimate of $p(\theta|y)$ will in some sense be less accurate than if you had integrated out $\alpha$ by hand and run the sampler on the marginal distribution. So the quality of the sampler's output isn't exactly the same, except at the limit. $\endgroup$
    – jbowman
    Jul 21 '18 at 14:47
  • $\begingroup$ @Xi'an - Well I did not know that, for some reason I had always believed it did! $\endgroup$
    – jbowman
    Jul 21 '18 at 14:49
  • $\begingroup$ @Xi'an - can you point me to a source for that? I simply cannot see how that could be the case; it seems to me you must lose information by generating a sample of $\theta$s conditional upon a sample of $\alpha$s relative to generating directly from the marginal distribution, even given that the $\alpha$s are generated from the conditional distribution. You would certainly lose information if the $\alpha$s were generated from their marginal distribution, after all... $\endgroup$
    – jbowman
    Jul 21 '18 at 15:01
  • $\begingroup$ In the marginal case, you can generate i.i.d. samples of the $\theta$, but in the Gibbs sampler case, they aren't i.i.d. any more, so for any given sample size, the estimates won't be as accurate, or so it seems to me. $\endgroup$
    – jbowman
    Jul 21 '18 at 15:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.