Hello Everyone,

I am implementing a marginalized poisson-gamma hierarchical model in pystan and I never get to initialize it properly. Any help will be appreciated.

Thanks

Model:

```
n ~ poisson(p)
c ~ gamma(n, m)
```

Marginalization

```
c ~ exp(-p) if c==0
c ~ sqrt(p/m*c) * exp(-c/m -p) * BessI_1(2*sqrt(c*p/m)) if c>0
```

Stan Code:

model_1 = “”“

functions {

real p_lpdf(real y, real p, real m){

if (y == 0) return (-p);

else {

return (0.5 * log(p / (y * m)) - y/m - p + log(bessel_first_kind(1, 2 * sqrt(y * p / m))));}

}

}

data {

int<lower=0> N;

vector[N] y;

}

parameters {

real<lower=0.1> lam;

real<lower=0.1> alpha;

}

model {

for (n in 1:N){

target += (p_lpdf(y[n]| lam, alpha));

}

}

”""

pooled_data_dict = {‘N’: 1000, ‘y’: y}

pooled_fit = pystan.stan(model_code=model_1,data=pooled_data_dict, iter=1000, chains=1)