Optimize models with arrays as data?

[JasonPekos]

Some models provide data in NxM arrays, e.g. the IRT model has (presumably participant x question) response data:

data {
  int<lower=0> I;
  int<lower=0> J;
  array[I, J] int<lower=0, upper=1> y;
}

Currently this is implemented like this:

for i in 1:I
    y[i,:] ~ arraydist([BernoulliLogit(a[i] .* (theta[j] .- b[i])) for j in 1:J])
end

But I feel like this isn't optimal. Leaving this here to revisit later.

[JasonPekos]

See also: dogs.