# Using Monte Carlo to sample from marginal distribution

I am defining a model on a vector, $$T$$, of size $$n$$, wherein each element $$t_i \in T$$ is independent and either $$0$$ or $$1$$. This model depends on 3 other parameters, $$q$$ (also a vector of size $$n$$), $$\theta$$, and $$\phi$$ (both $$\theta$$ and $$\phi$$ are vectors of size $$m$$). All elements of all vectors are independent. The vectors $$T$$ and $$q$$ have all been observed so I am going to assume that $$P(T)=c_1$$ and $$P(q) = c_2$$ where $$c_1$$ and $$c_2$$ are constants. Additionally, I assume a uniform prior on $$\theta$$, and some arbitrary prior on $$\phi$$ that is easy to sample from.

I define my model as follows:

\begin{align*} P\left(T\vert q,\theta,\phi\right)=\prod\limits_{k=1}^m\prod\limits_{\{i\mid h(i)=k\}}\left(q_i\theta_k\right)^{t_i}\left(1-q_i\theta_k\right)^{1-t_i}\tag{1} \end{align*} Where $$h$$ is some function that maps an index in $$T$$ to an index in $$\phi$$.

Since the $$t_i$$ are $$0$$ or $$1$$, this is essentially a product of bernoulli distributions each with a different parameter $$p=q_i\theta_k$$.

Ultimately, my goal is to estimate \begin{align*} P\left(\phi\vert T, q\right) \propto P\left(T\vert q,\phi\right)P\left(\phi\right)\tag{2} \end{align*}

The actual values of the $$\theta$$ are not relevant to me so I would like to take the marginal distribution. \begin{align*} P\left(T\vert q, \phi\right)=\prod\limits_{k=1}^m\int_1^0\prod\limits_{\{i\mid h(i)=k\}}\left(q_i\theta_k\right)^{t_i}\left(1-q_i\theta_k\right)^{1-t_i}d\theta_k\tag{3} \end{align*} This integral is quite difficult to do analytically (unless I can justify switching the order of the product and the integral) so my idea would be to try and sample directly from the distribution.

I am already doing MCMC to estimate the distribution in equation $$(2)$$ So I'd like to sample \begin{align*} \theta_k\sim U(0,1) \end{align*} Then plug this into equation $$(1)$$ and use that in my MCMC simulation for estimating equation $$(2)$$.

Does this sound like it would work? Is there a better way to estimate equation $${3}$$?

• There is no $\phi$ in (1)... Nov 12 '21 at 20:43
• @Xi'an, oops. thank you for pointing that out. The $\phi$ is encoded in the limit of the left most product and also in the function $h$ that maps indices of $T$ to indices of $\phi$. I'm not sure how best to notate this relation. Nov 13 '21 at 22:25

Using an unbiased estimate $$\hat\pi(\phi)$$ of the target density $$\pi(\phi)$$ in a Metropolis-Hastings algorithm is called the pseudo-marginal Metropolis–Hastings algorithm, first introduced by Beaumont (Genetics, 2003) and later extended by Andrieu and Roberts (Annals of Statistics, 2009). It is validated by treating $$\hat\pi(\phi)$$ as an auxiliary variable.