Hi everyone,

I need some help including discrete parameters in the model. I found some posts on Compound Poisson Gamma models but I couldn’t really follow. How can we implement this model in Stan?

In my case, the scale parameter of the Gamma distribution is N_j, which is Poisson distributed itself. I believe that I can’t write it as I am doing cause HMC is not optimized for discrete parameters. Since N_j is unbounded I couldn’t find a solution in the manual.

Is it possible to implement a continuous Poisson distribution?

f(x,\lambda)=a(\lambda)*\frac{e^{-\lambda}\lambda^x}{\Gamma(x+1)}

Here is my code:

```
data {
int<lower=0> N; //number of observations
vector<lower=0>[N] S_j; //Aggregate loss by j (response variable)
vector<lower=0>[N] N_j; //Aggregate number of observations by j
int<lower=0> K; //number of states
int<lower=0> M; //number of states-1
matrix[N,K] design_matrix_1;
matrix[K,M] design_matrix_2;
}
transformed data {
}
parameters {
// vector<lower=0>[K] lambda_j;
vector<lower=0>[K] theta_j;
real beta_0;
vector[M] beta;
vector<lower=0>[M] alpha_j;
real<lower=0> alpha_0;
real<lower=0> sigma_mu;
real<lower=0> sigma_alpha;
}
transformed parameters {
vector[N] theta_index = design_matrix_1 * theta_j;
// vector[N] lambda_index = design_matrix_1 * lambda_j;
vector[K] mu_index = beta_0 + design_matrix_2 * beta * sigma_mu;
vector[K] alpha_index = alpha_0 + design_matrix_2 * alpha_j * sigma_alpha;
}
model {
target += cauchy_lpdf(alpha_0 | 0, 5);
target += cauchy_lpdf(sigma_mu | 0, 1);
target += cauchy_lpdf(sigma_alpha | 0, 25);
target += normal_lpdf(alpha_j | 0, 1);
target += normal_lpdf(beta | 0, 1);
target += normal_lpdf(beta_0 | -9, 1);
target += gamma_lpdf(theta_j | alpha_index, alpha_index .* exp(-mu_index));
//target += poisson_lpmf(N_j | lambda_index);
target += gamma_lpdf(S_j | N_j , theta_index);
}
generated quantities{
real yrep[N]= gamma_rng( N_j , theta_index);
}
```

instead of