I’m trying to figure out a clever way to deal with a nasty integral inside my posterior.

I’m building a population model which has a Poisson-like term in the posterior to account for the estimated number of objects that should appear in a survey.

i.e.

Pr(\psi,\phi_{1:N}, z_{1:N} \vert N, d_{1:N}, \text{obs}, I) \propto Pr(\psi \vert I) e^{-\bar{N}(\psi, \text{obs.})}\prod_{i=1}^{N} n(z_i, \psi) Pr(\phi \vert z, \psi) C^{\text{DV}}(z_i) Pr(d_i \vert\phi_i, z_i)

Here, \psi are hyper parameters, \phi are object parameters and

\bar{N}(\psi, \text{obs.}) = \int d\Omega_{\phi} \int dz\; n(z_i, \psi) Pr(\phi \vert z, \psi) C^{\text{DV}}(z_i) Pr(S \vert \phi, z, \text{obs}).

The details aren’t so important, but I could not figure out a nice way to do this integral (~ 6 dimensions) numerically inside Stan. I thought that perhaps I precompute the integral outside Stan and then do a ~6D interpolation inside Stan to estimate \bar{N} for each draw on the hyper parameters, but interpolation seems to be tricky in Stan. Additionally, I can imagine all kinds of problems with going outside of the interpolation boundaries because I define my \psi as being normally distributed.

Unfortunately, these integrals are far from analytic thanks to physics and taking care of selections effects (Pr(S \vert \phi, z, \text{obs})).

Anyone run across some clever ways to handle this type of problem before?

Thanks!