Help implementing Marginalized Distribution

Modeling
1

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)

2

I don’t see, how they are connected, if you say n ~ poisson( c ), then n is distr. negative binomial.

3

I have just modified the post to represent what I meant. The actual model is:

Model:

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

Thanks for pointing out to my mistake.

4

Could you provide the derivation of the math behind? Some paper?

5

What is the story behind having a Poisson produce an integer parameter for a gamma?

Your model could run into problems if n == 0.

You can’t implement your model in Stan unless you can marginalize out n, i.e.,

p(c | m, p) = SUM_n poisson(n | p) * gamma(c | n, m)

Sometimes you can get away with restricting n to a range that’s tractable if p is small and all the probability mass lies near p.

Alternatively, you can use a continuous approximation

n ~ ... some positive continuous distribution
        with mean p and variance p ...