Hi there, I am a total newbie to Stan and I was wondering if anybody had any ideas on how to parallelise this code by David Phillipo to run over multiple cores? The issue I have is that the models run extremely slowly and often we have to batch process 50+ of these which soon becomes unmanageable. I wish I could share an attempt but I’m not even really sure where to begin other than some sense that I should be using reduce_sum(). There is a link to the code below and I’ve pasted the section of the code that I **think** can be parallelised underneath. Any help gratefully received, even if it’s just a hint to get started.

Link to code: multinma/binomial_2par.stan at master · dmphillippo/multinma · GitHub

Code section that I **think** can be parallelised:

```
// -- IPD likelihood --
if (link == 1) { // logit link
// Could replace with bernoulli_logit_glm in Stan > 2.20
ipd_r ~ bernoulli_logit(eta_ipd);
} else {
ipd_r ~ bernoulli(theta_ipd);
}
// -- AgD likelihood (arm-based) --
// We have to hand code the log likelihood contribution for the adjusted
// binomial here, as N is not necessarily an integer (which Stan doesn't
// like). The following is exactly equivalent to:
// ag_r ~ binomial(nprime, pprime);
for (i in 1:ni_agd_arm)
target += lchoose(nprime[i], agd_arm_r[i]) +
lmultiply(agd_arm_r[i], pprime[i]) +
(nprime[i] - agd_arm_r[i]) * log1m(pprime[i]);
}
generated quantities {
vector[ni_agd_arm * n_int_thin] theta2_bar_cum;
#include /include/generated_quantities_theta_fitted.stan
#include /include/generated_quantities_common.stan
#include /include/generated_quantities_theta.stan
// IPD log likelihood and residual deviance
for (i in 1:ni_ipd) {
log_lik[i] = bernoulli_lpmf(ipd_r[i] | theta_ipd[i]);
resdev[i] = -2 * log_lik[i];
fitted_ipd[i] = theta_ipd[i];
}
```