I am attempting to optimise a model which uses spherical harmonics to parameterise the success probability of binomials on the sphere. Transforming from the alm parameters to the logit-probability x involves two steps which are computationally expensive (the bracketed lines in the transformed parameters block). This operation will be the building block of a much more complicated model, so I’d like to check that I haven’t missed a trick. Can anyone suggest a different way of calculating this expression which is more Stan-friendly?

```
data {
int<lower=0> J; // number of pixels
int<lower=0> N; // number of HEALPix isolatitude rings
int<lower=0> H; // number of harmonics
int<lower=0> L; // 2 * max l of hamonics + 1
matrix[N,H] lambda; // spherical harmonics decomposed
matrix[J,L] azimuth; // spherical harmonics decomposed
int indices[J]; // map N->J
int lower[L]; // map H->L
int upper[L]; // map H->L
int k[J];
int n[J];
}
parameters {
vector[H] alm;
}
transformed parameters {
vector[J] x;
{
matrix[N,L] F;
for (i in 1:L) {
F[:,i] = lambda[:,lower[i]:upper[i]] * alm[lower[i]:upper[i]];
}
x = rows_dot_product(F[indices],azimuth);
}
}
model {
alm ~ std_normal();
k ~ binomial_logit(n, x);
}
```