Hi! I’m a fairly new Stan user and this is my first time posting to the forum so I hope I’m formulating everything correctly and apologies if not (also thank you to everyone who asks and answers questions on this group, it’s an amazing resource!).

I’m using @bgoodri’s excellent parameterisation of the multivariate probit model (original available here).

To make the program integrate more easily with `library(loo)`

and the new `reduce_sum`

parallelisation I tried to redefine the bulk of the model inside a single _lpmf function but I can’t work out how to vectorise the function to improve the computational performance.

The function as it stands is:

```
functions {
real multiprobit_lpmf(
int[] y, // response data - integer array of 1s and 0s
vector x, // covariate data
matrix beta, // coefficients for covariates
real[] u, // nuisance that absorbs inequality constraints
int D, // response dimension size
matrix L_Omega // lower cholesky of correlation matrix
){
vector[D] mu; // means
vector[D] z; // latent continuous state
vector[D] logLik; // log likelihood
real prev;
mu = beta * x;
prev = 0;
for (d in 1:D) { // Phi and inv_Phi may overflow and / or be numerically inaccurate
real bound; // threshold at which utility = 0
bound = Phi( -(mu[d] + prev) / L_Omega[d,d] );
if (y[d] == 1) {
real t;
t = bound + (1 - bound) * u[d];
z[d] = inv_Phi(t); // implies utility is positive
logLik[d] = log1m(bound); // Jacobian adjustment
}
else {
real t;
t = bound * u[d];
z[d] = inv_Phi(t); // implies utility is negative
logLik[d]= log(bound); // Jacobian adjustment
}
if (d < D) prev = L_Omega[d+1,1:d] * head(z, d);
// Jacobian adjustments imply z is truncated standard normal
// thus utility --- mu + L_Omega * z --- is truncated multivariate normal
}
return sum(logLik);
}
}
```

And the data being used:

```
data {
int<lower=1> K; // number of covariates + intercept
int<lower=1> D; // response dimension size
int<lower=0> N; // sample size
int<lower=0,upper=1> y[N,D]; // response data - integer array of 1s and 0s
vector[K] x[N]; // covariate data
}
```

At the moment the model is retrieving simulation parameters well but it would be nice to be able to remove the for loop, particularly when wrapping the function into partial_sums functions

```
model {
L_Omega ~ lkj_corr_cholesky(1);
to_vector(beta) ~ normal(0, 10);
// implicit: u is iid standard uniform a priori
{ // likelihood
for(n in 1:N){
target += multiprobit_lpmf(y[n] | x[n], beta, u[n], D, L_Omega);
}
}
}
```

Any advice on how to vectorise the function (or if this is possible/desirable) is greatly appreciated- thanks!