# MCMC when the density involves integration over a simplex

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$.]

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,