I have a vector of observations y (of length N) with associated indices (ii, jj, kk), each of length N, and an array of parameters theta[i, j, k] in my model, and I would like to write a vectorised sampling statement equivalent to
for (l in 1:N) y[l] ~ normal(theta[ii[l], jj[l], kk[l]], sigma);
Is it possible to vectorise this directly? If not, does it make sense to do it indirectly, somehow, like vectorizing (“raveling”) the three-dimensional array and computing the indices to the vector manually?
The only comment I would make is that, if this is the only time theta is used, then could theta be a single dimensional vector of length L? In that case vectorization would be possible. Is theta being used elsewhere?
For reference (if it helps), 3 dimensional arrays are C++ std::vectors of std::vectors of std::vectors. Stan Matrix and vector types are Eigen Matrix and Vector types. The biggest performance boost from vectorization is when there’s room for lots of shared calculation. I don’t think that’s true for a bunch of scalar normal_lpdf evaluations, but I don’t have a great feel for what’s gonna be fast and slow here. You could very well get some performance boost by shuffling things around and making the memory happier.