Yeah, the code could be made more efficient. One more sketch to illustrate some compression (I hope all of those work):

```
transformed data {
matrix[N, num_options] skipped_options = rep_matrix(0, N, num_options);
for(n in 1: N) {
for(option in 1:(observed_option[n] - 1)) {
skipped_options[n,option] = 1;
}
}
}
model {
vector[N] val = beta * X;
matrix[N,num_options] gamma_rep = skipped_options .* rep_matrix(to_rowvector(gamma), N);
matrix[N,num_options] val_rep = skipped_options .* rep_matrix(val, num_options);
matrix[N,num_options] skipped_options_lp = skipped_options .* log1m(normal_cdf(gamma_rep - val_rep, 0, 1));
target += to_vector(skipped_options_lp);
target += normal_lcdf(gamma[observed_option[n]] - val, 0, 1);
}
```

I guess I see what you meant by the “latent binaries” (the `skipped_options`

matrix) but those can be computed in transformed data so should not make problems. There will be some wasted computation in the matrix operations so you should test whether this is an improvement, but I guess it could be worth it. I also multipled `gamma_rep`

and `val_rep`

with `skipped_options`

in hope it reduces the autodiff stack (the elements multiplied by zero won’t propagate) but this is also just a guess and might not be an improvement in practice.

Another approach would be to group the data by the value of `observed_option`

which will leave you with an iteration over options and the matrices involved at each step won’t have any zeroes and no computation will be wasted.