MCMC combined with numerical integration towards more efficient Bayesian inference I am quite new to Bayesian statistics so the question can be a bit naive.
My question is on how to deal with a model with individual coefficients. Simple versions of a task and a model I deal with is bellow.
Let there are customers that visit our company from time to time. We are interested in modeling the inter-visit time intervals.
The intervals follow the gamma distribution:
$$
t_{ij} \sim \mathtt{gamma}(k,\theta),
$$
where $i$ is the customer number and $j$ the time interval number for the customer $i$ with $k \sim N(0,10)$.
However we also assume that customers are heterogeneous, which means that every customer has its own individual parameter $\theta_i$ and, in fact, we have the following model:
$$
t_{ij} \sim \mathtt{gamma}(k,\theta_i),
$$
where $\theta_i$ comes from the inverse gamma distribution:
$$
\frac{1}{\theta_i} \sim \mathtt{gamma}(\nu,\sigma),
$$ with $\nu,\sigma \sim N(0,10)$. Thus, the model has 3 parameters: $k,\nu, \sigma$.
The question is how to encode this model (in particular in stan) for efficient inference? The number of observations is high, e.g. 10000+.
The first option is to encode it straightforward and allow MCMC sample a large space including $\theta_i$ for every customer. However, it is at least slow and looks like the number of iteration in MCMC should be high.
Would it be correct to remove $\theta_i$ from the sampling space and substitute it with numerical integration for computing the target log likelihood as follows:
$$
t_{ij} \sim \int \mathtt{gamma}(k,\theta)\mathtt{gamma}\left(\frac{1}{\theta} | \nu,\sigma\right)d{\theta},
$$
where the integral can be computed by taken for example just 1000 (instead of 10000+) points of $\theta^{(i)}$ (given $\nu$ and $\sigma$) that are not included into the MCMC sampling space.
PS: If "model with individual coefficients" is an incorrect term, I would be grateful to know the correct one, since the correct terms dramatically increase the search possibilities.
UPD: I describe the exact model I want to implement, since without it my question seems to be misguiding. The model is from Allenby G.M., Leone R.P., Jen L. A Dynamic Model of Purchase Timing with Application to Direct Marketing // J. Am. Stat. Assoc. 1999. Vol. 94, № 446. P. 365. For simplicity I provide no priors.
$$
t_{ij} \sim \Phi(\beta_i\cdot x_{ij})\cdot \mathtt{ggamma}(k_1,\theta^{(1)}_i,\gamma) + (1-\Phi(\beta_i\cdot x_{ij}))\cdot \mathtt{ggamma}(k_2,\theta^{(2)}_i,\gamma),
$$
$$
\frac{1}{\theta^{(1)}_i} \sim \mathtt{ggamma}(\nu_1,\psi_1,\gamma),
$$
$$
\frac{1}{\theta^{(2)}_i} \sim \mathtt{ggamma}(\nu_2,\psi_2,\gamma),
$$
$$
\beta \sim N(\mu,\sigma),
$$
where $\Phi(x)$ is normal CDF. So the parameters of the model are $\mu,\sigma,$ (vectors corresponding to the covariates $x_{ij}$) $\nu_{1,2},\psi_{1,2}, k_{1,2}$ (just 6 numbers). The value of $\gamma$ is considered as a constant and is not from the optimization space.
 A: A first reduction of complexity is to use sufficiency: if
$$t_{ij} \sim \mathtt{gamma}(k,\theta_i),$$
then
$$t_{i\cdot} = \sum_{j=1}^J t_{ij} \sim \mathtt{gamma}(J\times k,\theta_i),$$
is sufficient. And integrating out the $\theta_i$'s is also feasible by conjugacy
$$\int_0^\infty \theta_i^{-Jk-\nu-1}\exp\{-\theta_i^{-1} [t_{i\cdot}+\sigma]\}\,\text{d}\theta_i = \Gamma(Jk+\nu)[t_{i\cdot}+\sigma]^{Jk+\nu} $$
So there is no need for a call to simulation methods.
A: What you have in the OP is not exactly true, but there is a correct way to do this sort of thing. Suppose that I have parameters $\theta_1, \ldots, \theta_N \sim g(\theta \mid\eta)$ given $\eta$, and I use a prior $\eta \sim F$. Then the joint distribution of the $\theta$'s after integrating out $\eta$ is 
$$
m(\pmb \theta) = \int \left(\prod_{i = 1}^N g(\theta_i \mid \eta) \right) \, F(d\eta). \tag{$\dagger$}
$$
A "folk-theorem" in MCMC is that chains based on $(\dagger)$ leads to better mixing, i.e., it is better to marginalize things out if you can. It is not strictly true, as there are some situations where parameter-expansion ideas lead to better mixing, and it also matters what exactly you do with $m(\pmb \theta)$ to sample it, but it is a useful rule of thumb. 
So, for example, it would be valid, in STAN to do 
target += log_m(theta)

assuming you can write a function log_mthat computes the log of $(\dagger)$. That's the tricky part, of course, since log_m may not be easy to evaluate to a given accuracy. 
Actually, this is how the makers of STAN recommend that one incorporates discrete parameters, because HMC cannot deal with discrete parameters directly. They recommend summing out the discrete parameter, which is equivalent to $(\dagger)$ when $\eta$ is discrete (e.g., it might be a mixture component indicator). 
