I’m still struggling here. Do I really want to increase target with both pi as well as with zeta? For some reason that seems off to me, but I’m pretty new at this so it’s quite possible.

When I use (a multivariate version of) the model below I get divergent sample errors. Is this model that you’re proposing? If so, I’ll look into issues with incorporating my additional feature variables, else I’ll try to fix my model.

```
data {
int<lower=1> M; // num observations
int<lower=1> N; // number of sub-observations for all observations
real<lower=0> y[N]; //Observations where each row is a subobservation
int<lower=1,upper=M> observation[N]; //observation id for subobservation n
int<lower=1> row_to_subobservation_map[N]; //subobservation index for each row
int<lower=1> count_of_subobservations_per_observation[M];
real shape; // subobservation weight prior
}
transformed data {
vector<lower=0>[N] shape_vec;
for (n in 1:N){
shape_vec[n] = shape; //symmetric dirichlet prior
}
}
parameters {
vector<lower=0>[N] zeta; // subobservation dist over observation
}
transformed parameters {
}
model {
for (n in 1:N) {
int zeta_mark = n-row_to_subobservation_map[n]+1;
int nc = count_of_subobservations_per_observation[observation[n]];
vector[nc] pi = segment(zeta,zeta_mark,nc); // gamma(zeta | shape, 1)
pi = pi / sum(pi); // thus dirichlet(pi | shape)
target += log(y[n]*pi[row_to_subobservation_map[n]]); // likelihood
}
target += gamma_lpdf(zeta | shape, 1.);
}
```