I am interested in samples of $\theta$ from the posterior distribution

$$ P(\theta|x) = \int d\phi P(\theta|\phi)P(\phi|x) $$

where $x$ are data and $\phi$ are nuisance parameters. In principle, I can use a Metropolis Hastings sampler to sample from $\theta$ and $\phi$ and discard the samples of the nuisance parameter.

In this case, I can sample from $P(\phi|x)$ directly such that I can approximate the marginal posterior by Monte Carlo integration

$$ P(\theta|x)\approx\frac{1}{n}\sum_{i=1}^n P(\theta|\phi_i)\equiv \hat P, $$

where $\phi_i$ are samples from $P(\phi|x)$. Of course, the approximation $\hat P$ is not deterministic because of the sampling error. I seem to remember that running the following algorithm samples from the posterior in this case but I can no longer find the reference. Is this correct? If so, do you have a reference?

# Sampler in pseudo-python
theta = initial_value

for step in range(num_steps):
    # Sample phi
    phi = sample_phi_given_x(x)
    # Evaluate current posterior
    posterior_estimate = mean(posterior(theta, phi))
    # Sampled from proposal and evaluate posterior at proposal
    candidate = sample_from_proposal(theta)
    posterior_candidate_estimate = mean(posterior(candidate, phi))

    # Accept or reject the proposal
    if random_uniform() < posterior_candidate_estimate / posterior_estimate:
         theta = candidate
  • $\begingroup$ This classic form of Markov Chain Monte Carlo produces only well-defined posterior distributions. $\endgroup$ Commented Mar 30, 2017 at 14:57
  • $\begingroup$ @MichaelChernick, I'm afraid I don't follow. $\endgroup$ Commented Mar 30, 2017 at 14:59
  • 3
    $\begingroup$ Look up here: darrenjw.wordpress.com/2010/09/20/… (for a reference, the term you need to look for is pseudo-marginal MCMC). $\endgroup$
    – lacerbi
    Commented Mar 30, 2017 at 15:19
  • $\begingroup$ @lacerbi: thanks for the link! @eric_kernfeld: an independence sampler is not suitable in my case because I don't have a good idea of the high-density regions of the posterior. Removing the call to mean will reduce the acceptance rate if the variance of the posterior associated with sampling $\phi$ is sufficiently large. In other words, the algorithm is hit by the curse of dimensionality. $\endgroup$ Commented Mar 30, 2017 at 15:32

1 Answer 1


I do not understand the need for this construction if you can

  1. simulate from $P(\phi|x)$
  2. simulate from $P(\theta|\phi)$

since neither an approximation nor an MCMC implementation is then required.

If you cannot analytically sample from P(θ|ϕ), a regular Gibbs sampler will work as well:

  1. simulate from $P(\phi|x,\theta)\propto P(\phi|x)\times P(\theta|\phi)$ [which may require one [& only one] Metropolis-within-Gibbs step]
  2. simulate from $P(\theta|\phi)$ [which may require one [& only one] Metropolis-within-Gibbs step]
  3. simulate from $P(\phi|x,\theta)\propto P(\phi|x)\times P(\theta|\phi)$ &tc...

This is less costly than the pseudo-marginal version, which requires $n$ $\phi_i$'s for one proposal of $\theta$.

Note: The link to pseudo-marginal MCMC made by lacerbi is however correct in that if one manages to get an unbiased estimator of the posterior density (up to a constant) there are valid MCMC implementations based on such estimates. Which means your proposal is not correct because you use the same $\phi_i$'s for numerator and denominator in the Metropolis ratio. The random variables used to approximate the posterior density should be considered as auxiliary or latent variables and hence be only simulated for the proposal, while being preserved from the earlier stage for the current value: to quote from Darren Wilkinson's blog,

"The key to understanding the pseudo-marginal approach is to realise that at each iteration of the MCMC algorithm a new value of w is being proposed in addition to a new value for x."

Meaning that the acceptance ratio

enter image description here

keeps the previous value of w for the previous value of x. In other words, the $\phi_i$'s are not to be re-simulated for computing the unbiased estimator at the previous value of $\theta$.

  • $\begingroup$ Good point. In practice, I'd like to avoid the two-stage sampling because $P(\theta|\phi)$ has high variance with respect to $\phi$ and I cannot analytically sample from $P(\theta|\phi)$. So I'd have to (1) draw a sample of $\phi$, (2) run an MH chain with target distribution $P(\theta|\phi)$ to convergence, and repeat. Or is there a better approach? $\endgroup$ Commented Mar 31, 2017 at 12:33
  • $\begingroup$ In the proposed solution, $\hat P = \frac 1n \sum_{i=1}^n P(\theta|\phi_i)$ has variance wrt $\phi$ which is smaller than the variance of $P(\theta|\phi)$ by a factor of $n$. $\endgroup$ Commented Mar 31, 2017 at 12:44
  • $\begingroup$ I think the algorithm I proposed above is correct because the distribution is reevaluated on the same sample for the current state and the candidate state. Details here. $\endgroup$ Commented Mar 31, 2017 at 15:02
  • $\begingroup$ Thank you for all the insights. The problem with my sampling was a result of the estimate of the posterior not being unbiased. I will probably go for the Gibbs sampling approach. $\endgroup$ Commented Apr 1, 2017 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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