Vectorize multinomial_lpmf

Currently, there is only this signature:
real multinomial_lpmf(array[] int y | vector theta) (the vector is actually a simplex but this is probably for efficiency?)

I wish there was a vectorized variant with the following signature (or compatible):

real multinomial_lpmf(array[ , ] int y | array[] simplex theta)

and probably

real multinomial_lpmf(array[ , ] int y | simplex theta) (assume all have same parameters)

it feels awkward to write for cycle for something most other distributions implement.

sample code where I use it to estimate mean product rating:

data {
  int<lower=1> N;  // num products
  array[N, 5] int<lower=0> ratings;
}
parameters {
  vector<lower=0>[5] alpha;
  array[N] simplex[5] d;
}
model {
  d ~ dirichlet(alpha);  // for all d[i]
  for (i in 1:N) {  // unlike dirichlet, multinomial not vectorized
    ratings[i] ~ multinomial(d[i]);
  }
}

Start at the multinomial_lpmf entry point and compare its current signature with vectorized distribution interfaces in the library. Define compatible overload behavior for arrays of counts and simplexes, including the shared-parameter case, and verify the sample model no longer needs an explicit loop.