# Help implementing Marginalized Distribution

**URL:** <https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550>\
**Category:** Modeling\
**Created:** [August 12, 2017, 10:31pm UTC](https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550 "2017-08-12T22:31:45Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![mig](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mig/32/14875_2.png) [@mig](https://discourse.mc-stan.org/u/mig)\
**Post date:** [August 12, 2017, 10:31pm UTC](https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550/1 "2017-08-12T22:31:45Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [August 13, 2017, 1:46am UTC](https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550/2 "2017-08-13T01:46:41Z")

</div>

> [@mig](#):
>
> Model:
> 
> ```
> n ~ poisson(p)
> c ~ gamma(p, m)
> 
> ```

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

---

<div class="post-metadata">

**Author:** ![mig](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mig/32/14875_2.png) [@mig](https://discourse.mc-stan.org/u/mig)\
**Post date:** [August 13, 2017, 1:45pm UTC](https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550/3 "2017-08-13T13:45:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [August 14, 2017, 4:14am UTC](https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550/4 "2017-08-14T04:14:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [August 15, 2017, 7:48pm UTC](https://discourse.mc-stan.org/t/help-implementing-marginalized-distribution/1550/5 "2017-08-15T19:48:32Z")

</div>

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 ...

```