Help Deriving Likelihood Term When the Target is Known Probabilistically I am trying to model data $\{Y_t,Q_t\}_{t=1}^T$, where the model is parameterized by $\theta$. $Y_t$ is a quantity where the model prediction can be solved in closed form, $\hat{Y}_t(\theta)$, where the model prediction of $Q_t$, $\hat{Q}_t(Y_t,\theta)$, can only be simulated via Monte Carlo. The simulation results in an estimate of the mean, $\bar{Q}_t(Y_t,\theta)$, and Monte Carlo variance, $\sigma_t(Y_t,\theta)^2$. Thus, $\hat{Q}_t(Y_t,\theta) \sim N\left(\bar{Q}_t\left(Y_t,\theta\right),\sigma_t\left(Y_t,\theta\right)^2\right)$.
I would like to compute the likelihood. For a specific $t$, the contribution to the likelihood is $$p(Y_t, Q_t | \theta) = p(Y_t|\theta)p(Q_t|Y_t, \theta)$$.
How do I incorporate the uncertainty of $\hat{Q}_t(Y_t,\theta)$ caused by the Monte Carlo simulation in the likelihood? Setting $\hat{Q}_t(Y_t,\theta) = \bar{Q}_t(Y_t,\theta)$ would ignore the uncertainty.
My instinct is to try to integrate out the noise, but I am not sure if that is technically correct:
$$p(Y_t, Q_t | \theta) = p(Y_t|\theta)p(Q_t|Y_t,\theta) = p(Y_t|\theta) \int_{\mathbb{R}}p\left(Q_t| \hat{Q}_t\left(Y_t, \theta\right) = X,Y_t, \theta\right)p\left(\hat{Q}_t\left(Y_t,\theta\right) = X | Y_t,\theta\right) dX$$

I believe I was correct in that assuming that I can marginalize out the noise from the joint-density of $Q_t$ and $\hat{Q}_t$
$$p\left(Q_t \middle| Y_t, \theta\right) = \int_{\mathbb{R}} p\left(Q_t, \hat{Q}_t\left(Y_t,\theta\right) = X \middle| Y_t,\theta\right)dX$$.
While it wasn't explicitly stated previously,
$$p\left(Q_t, \hat{Q}_t\left(Y_t,\theta\right) \middle| Y_t,\theta\right) = p\left(Q_t \middle| \hat{Q}_t\left(Y_t, \theta\right), Y_t , \theta\right) \times p\left(\hat{Q}_t\left(Y_t,\theta\right) \middle| Y_t, \theta\right)$$
follows from Bayes' Law.
It should be noted sometimes these integrals have closed-form solutions. For example when $p\left(Q_t \middle| \hat{Q}_t\left(Y_t, \theta\right), Y_t , \theta\right)$ is also Gaussian. Which is the case encountered many times where a Gaussian prior is assumed in Bayesian statistics. In my specific example the conditional density is log-Gaussian, numerical integration techniques such as Gaussian quadrature must be used.
