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}$?

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.
