# 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.

• I can't really help you with this problem, but for details on coding the model in Stan it might be faster to ask on the Stan forums (discourse.mc-stan.org). Sep 21, 2018 at 7:13
• Thank you, Maurits! I will duplicate my question also there Sep 21, 2018 at 7:19
• It is not correct to do that because the $t_{ij}$’s are no longer independent when $\theta_i$ is integrated out.
– guy
Sep 23, 2018 at 13:50
• @guy, intuitively it is not correct, however, if I am not wrong, formally the result is the same. You just compute the probability the get all several values together that is basically the same as get them independently. The dependence is only needed to determine the $\theta_i$. Sep 23, 2018 at 14:10
• @guy, see alse the answer by Xi'an Sep 23, 2018 at 14:23

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.

• Thank you for this answer. However, it seems, that this is only valid if we consider gamma distribution. In the actual model there is a mixture of two generalized gamma distribution with individual $\theta$ as well as with individual probability for the mixture that depends on some covariates. You also cannot do the first step unfortunately since the weights in the mixture depend on the interval. Here I have simplified the task, probably two much. My main concern was on what can we do the individual coefficients. Sep 23, 2018 at 14:17
• If you ask a question about one model and want an answer about another model, what is the point of asking the question? Sep 23, 2018 at 14:37
• Sorry for this, my point is not to be too specific. I want to deal with individual coefficients, and the model is given as an example. If I would have asked about my model, it is unlikely that it would be useful for a large number of people, is not it? I am new to the community, so if it is better to provide the most specific description of the model, even if it is useful for only a few people, I will do that. Sep 24, 2018 at 4:28

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).

• Thank you! So it means that I can try to numerically integrate it. The most "fast-and-dirty" approach is a random sampling, i.e., Monte-Carlo method. Is it possible to randomly sample from a distribution in the model (I guess I should not put it in the parameters, since I do not want to optimize w.r.t. these sampling)? As far as I know *_rng functions are not available in the model part. Sep 24, 2018 at 4:50
• @Aleksey there are enough pitfalls in doing this type of thing numerically that I’m not sure I would recommend it to a novice. Monte Carlo integration will induce bias, and correcting this bias is non-trivial and probably cannot be done using STAN. And numerical integration is also not so easy in that, if $N$ is large, the mass is likely to be concentrated in a very tiny area (it is also quite expensive potentially).
– guy
Sep 24, 2018 at 15:35