Ragged array of simplexes

Modeling
4

That is quite possibly correct but more complicated than necessary. Dealing with simplexes of different sizes is annoying, but it it actually the least annoying of the ragged arrays. The trick to avoiding having to deal with Jacobians is to remember that a Dirichlet distribution over simplexes can be constructed by normalizing a set of unit-scale but independent Gamma random variables:


Thus, in Stan you declare a looooong vector of non-negative parameters

parameters {
  vector<lower=0>[sum(K)] gamma; // K is an int[] declared in data
  ...
}

and in the model block you break off however many of those you need at the moment, divide them by their sum, and use the resulting simplex in your likelihood, as in

model {
  int pos = 1;
  for (j in 1:J) {
    vector[K[j]] pi = segment(gamma, pos, K[j]);
    pi = pi / sum(pi); // simplex
    target += // some function of pi
    pos = pos + K[j];
  }
}

Then you put a Gamma prior on gamma

  target += gamma_lpdf(gamma | alpha, 1);
   // implies pi ~ dirichlet(segment(alpha, pos, K[j]))

where alpha is a looooong vector of known shape hyperparameters for the Dirichlet distributions that are concatenated together. If alpha = rep_vector(1, sum(K)); then each pi is uniform a priori over simplexes of size K[j] and you can get away with

  target += exponential_lpdf(gamma | 1);
  // implies pi ~ uniform over simplexes
3 Likes