Hello. I implemented a compound Generalized Dirichlet Multinomial, and the correponding `_rng`

function.

Question: When I sample from the model, and use that data as training data, in an attempt to infer the parameters, I cannot seem to recover them. Moreover, if I sample from the distribution using the returned parameters, the results are not similar to the training data.

Perhaps there is a bug in the model? I cannot seem to find literature to compute the `_rng`

so I used the generative process as above. Maybe that is the issue?

Any assistance would be greatly appreciated!

```
functions {
real dirichlet_multinomial_lpmf(int[] y, vector alpha) {
real alpha_plus = sum(alpha);
return lgamma(alpha_plus) + sum(lgamma(alpha + to_vector(y)))
- lgamma(alpha_plus+sum(y)) - sum(lgamma(alpha));
}
int[] dirichlet_multinomial_rng(vector alpha, int N) {
return multinomial_rng(dirichlet_rng(alpha), N);
}
real generalized_dirichlet_multinomial_lpmf(int[] y, vector alpha, vector beta_) {
// y is num_categories dimensional, alpha and beta are num_categories-1 dimensional
int D = dims(alpha)[1]; // D = num_categories-1
vector[D+1] x = cumulative_sum(to_vector(y));
vector[D] z;
vector[D] alpha_prime = alpha + to_vector(y)[1:D];
vector[D] beta_prime;
z[1] = x[D+1];
z[2:D] = rep_vector(x[D+1],D-1) - x[1:D-1];
beta_prime = beta_ + z;
return (lgamma(x[D]+1) - sum(lgamma(to_vector(y)+rep_vector(1,D+1)))
+ sum(lgamma(alpha_prime)) - sum(lgamma(alpha))
+ sum(lgamma(beta_prime)) - sum(lgamma(beta_))
- sum(lgamma(alpha_prime + beta_prime)) + sum(lgamma(alpha+beta_))
);
}
int[] generalized_dirichlet_multinomial_rng(vector alpha, vector beta_, int N) {
int D = dims(alpha)[1];
vector[2] tmp;
int out[D+1];
out[1] = N;
for(n in 1:D){
tmp[1] = alpha[n];
tmp[2] = beta_[n];
out[n:n+1] = dirichlet_multinomial_rng(tmp,out[n]);
}
return out;
}
}
data{
int<lower=1> num_obs; // number of observations
int<lower=1> num_test_obs; // number of test observations
int<lower=1> num_categories; // number of categories
int<lower=0> obs[num_obs,num_categories]; // data
int<lower=0, upper=1> run_estimation;
}
parameters {
vector<lower=0>[num_categories-1] alpha_std;
vector<lower=0>[num_categories-1] beta_std;
}
transformed parameters {
vector<lower=0>[num_categories-1] alpha;
vector<lower=0>[num_categories-1] beta_;
alpha = exp(1.0 + alpha_std);
beta_ = exp(1.0 + beta_std);
}
model{
alpha_std ~ normal(0,1);
beta_std ~ normal(0,1);
if (run_estimation==1){
for(n in 1:num_obs){
obs[n] ~ generalized_dirichlet_multinomial(alpha,beta_);
}
}
}
generated quantities {
int y_test[num_test_obs,num_categories]; // test data
if (run_estimation==0){
for (n in 1:num_test_obs) {
y_test[n] = generalized_dirichlet_multinomial_rng(alpha,beta_,1000);
}
}
}
```