I have the following setup.

  • Parameters $W$ with density $\pi(w)$.
  • Observed data $X_1,...,X_n$ iid.
  • Density of $X_i|W=w$ is $f(x_i|w) = \int_{\Delta(x_i)} f(\mathbf c|w) \,d\mathbf c$.
  • The simplex $\Delta(x_i) = \{\mathbf c \geq 0 : c_1 + \cdots + c_m = x_i\}$.

I want to find the distribution of $W|X_1=x_n,\dots,X_n=x_n$.

Although $f(\mathbf c|w)$ is a fairly cheap to evaluate, I don't think it is possible (after a lot of effort!) to get a closed or simple form for $f(x_i|w)$. Therefore it's also difficult to find out much about the posterior distribution of $W$, analytically, either.

Instead I have been considering MCMC for this purpose. It is straightforward (assuming I can find a good enough $q$) to write down a basic Metropolis-Hasting algorithm:

  • Come up with a proposal distribution $q(w|w')$ and pick $w_0$.
  • For $t=1,2,3,\dots$:
    • Sample $w$ from a $q(\cdot|w_{t-1})$.
    • Let $r = \frac{ q(w|w_{t-1}) \pi(w) \prod_{i=1}^n f(x_i|w) } { q(w_{t-1}|w) \pi(w_{t-1}) \prod_{i=1}^n f(x_i|w_{t-1}) }$.
    • Let $w_t = w$ with probability $\min(r,1)$, and let $w_t = w_{t-1}$ otherwise.

The trouble is that, as was already mentioned above, the only way I know to evaluate $f(x_i|w)$ is via a (multi-dimensional) numerical integration. And in the above algorithm, this has to be done $2n$ times per iteration!

So my $\textbf{questions}$ are, in order of specificity:

  • Is there a way to combine the numerical integration required for $f(x_i|w)$ and the overall MCMC algorithm? (Intuitively, suppose we use an MCMC algorithm to evaluate each integral $f(x_i|w)$, then there would be two levels of nested Markov chains, and perhaps there is a way to combine these chains to get better overall performance.)
  • Is there a (obvious, standard, or otherwise) better way to solve my problem than the basic MH algorithm above?

Thanks in advance for any replies!

[P.S. In case it matters, $f(\mathbf c|w)$ is the density of a probability distribution on $\mathbb R^m$.]


1 Answer 1


You could push the MCMC machinery one step further by simulating the $\mathbf{c}_i$'s conditional on $x_i$ and $w$. Indeed, your equation $$f(x_i|w)=\int_{\Delta(xi)}f(\mathbf{c}_i|w)d\mathbf{c}_i$$ states that $x_i$ is obtained by marginalisation over the $\mathbf{c}_i$'s. In statistical terms, the $\mathbf{c}_i$'s are latent variables: they are not observed, but they are useful in defining the sampling model.

Therefore, an alternative solution is to incorporate all the $\mathbf{c}_i$'s within your model and to simulate them as well as $w$ conditional on the $x_i$'s:

  1. Gibbs sampling version:

    • simulate the $\mathbf{c}_i$'s conditional on $x_i$ and $w$;
    • simulate $w$ conditional on the $\mathbf{c}_i$'s (since they are independent of the $x_i$'s, conditionally on the $\mathbf{c}_i$'s).
  2. Metropolis-Hastings version:

    • Make a proposal for a new value of $w$, e.g. $w\sim q(w|w_{t-1})$;
    • Given $w$ and the $x_i$'s, simulate the $\mathbf{c}_i$'s;
    • Accept the proposed $w$ with probability $$ \dfrac{q(w_{t-1}|w)}{q(w|w_{t-1})}\, \dfrac{\pi(w)\prod_i f(\mathbf{c}_i|w)}{\pi(w_{t-1})\prod_i f(\mathbf{c}_{i,t-1}|w_{t-1})}\,\wedge 1 $$

This way, you do not introduce an extra degree of approximation in the analysis. Even though simulating the $\mathbf{c}_i$'s obviously increases your computing time,


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.