# Zero Inflated Poisson Log-Lik

**URL:** <https://discourse.pymc.io/t/zero-inflated-poisson-log-lik/2664>\
**Category:** Questions\
**Created:** [February 6, 2019, 12:49am UTC](https://discourse.pymc.io/t/zero-inflated-poisson-log-lik/2664 "2019-02-06T00:49:44Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Brian](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/brian/32/7088_2.png) [@Brian](https://discourse.pymc.io/u/Brian)\
**Post date:** [February 6, 2019, 12:49am UTC](https://discourse.pymc.io/t/zero-inflated-poisson-log-lik/2664/1 "2019-02-06T00:49:44Z")

</div>

I am trying to learn Pymc3 by creating my own version of zip log-likelihood. I eventually want to expand to a case not covered by built in Pymc3.

If I run the same observed y’s through my implementation and the built in one, the sign of ‘ap’ is flipped (negative in the built-in, positive in mine). It seems like this should be straightforward to take the log of the likelihood:

f(x \mid \psi, \theta) = \left\{ \begin{array}{l} (1-\psi) + \psi e^{-\theta}, \text{if } x = 0 \\ \psi \frac{e^{-\theta}\theta^x}{x!}, \text{if } x=1,2,3,\ldots \end{array} \right.

Is there a simple bug in my code??

#MINE

```
class ZeroInflatedPoisson_JA(Continuous):
    def __init__ (self, pi, theta, *args, **kwargs):
        super(ZeroInflatedPoisson_JA, self). __init__ (*args, **kwargs)
       
        #create variables for theta and pi
        self.theta = theta = tt.as_tensor_variable(theta)
        self.pi = pi = tt.as_tensor_variable(pi)
        self.poisson_ = Poisson.dist(theta)

    def logp(self, value):
        return tt.switch(tt.eq(value, 0),
                         tt.log((1 - self.pi) + self.pi*tt.exp(-1.0*self.theta)),
                         tt.log(self.pi)+self.poisson_.logp(value)
                        )

with pm.Model() as ZIPP_JA:

    #prior
    ap = pm.Normal('ap', 0., 1.)
    al = pm.Normal('al', 0., 10.)
    
    
    #logit link (rethinking code is logit(p) <-ap which implies that p is invlogit(ap)
    #sigmoid is same thing as invlogit)
    p = pm.math.sigmoid(ap)
    theta = pm.math.exp(al)
    
    y_obs = ZeroInflatedPoisson_JA('y_obs', p, theta, observed=y)

    ZIPP_JA = pm.sample(1000, tune=1000,chains=4)
```

---

<div class="post-metadata">

**Author:** ![lucianopaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/lucianopaz/32/2186_2.png) [@lucianopaz](https://discourse.pymc.io/u/lucianopaz)\
**Post date:** [February 6, 2019, 11:06am UTC](https://discourse.pymc.io/t/zero-inflated-poisson-log-lik/2664/2 "2019-02-06T11:06:00Z")

</div>

To start, I don’t see any difference in `ap` between your implementation and pymc3’s current master branch. The only problems I do see are:

1. Your custom zero inflated poisson class should inherit from `Discrete` instead of from `Continuous`
2. You should check in your custom `logp` that `all(value >= 0)`, and that your parameter values are in the proper ranges (`[0, 1]` for `pi` and `[0, +inf)` for theta), and return `-inf` otherwise.

I can only recommend that you upgrade your pymc3 version.

> **This is the code I ran to compare both implementations**
>
> ```python
> import pymc3 as pm
> from pymc3.distributions import (Discrete, Poisson,
> ZeroInflatedPoisson)
> import theano.tensor as tt
> 
> class ZeroInflatedPoisson_JA(Discrete):
> def __init__ (self, pi, theta, *args, **kwargs):
> super(ZeroInflatedPoisson_JA, self). __init__ (*args, **kwargs)
> # create variables for theta and pi
> self.theta = theta = tt.as_tensor_variable(theta)
> self.pi = pi = tt.as_tensor_variable(pi)
> self.poisson_ = Poisson.dist(theta)
> 
> def logp(self, value):
> return tt.switch(tt.eq(value, 0),
> tt.log((1 - self.pi) +
> self.pi*tt.exp(-1.0*self.theta)),
> tt.log(self.pi)+self.poisson_.logp(value)
> )
> 
> P = 0.9
> THETA = 5
> y = ZeroInflatedPoisson.dist(psi=P, theta=THETA).random(size=1000)
> 
> with pm.Model() as ZIPP_JA:
> # prior
> ap = pm.Normal('ap', 0., 1.)
> al = pm.Normal('al', 0., 10.)
> p = pm.Deterministic('psi', pm.math.sigmoid(ap))
> theta = pm.Deterministic('theta', pm.math.exp(al))
> 
> y_obs = ZeroInflatedPoisson_JA('y_obs', p, theta, observed=y)
> 
> ZIPP_JA = pm.sample(1000, tune=1000, chains=4)
> pm.traceplot(ZIPP_JA)
> 
> with pm.Model() as pymc3_model:
> # prior
> ap = pm.Normal('ap', 0., 1.)
> al = pm.Normal('al', 0., 10.)
> p = pm.Deterministic('psi', pm.math.sigmoid(ap))
> theta = pm.Deterministic('theta', pm.math.exp(al))
> 
> y_obs = ZeroInflatedPoisson('y_obs', p, theta, observed=y)
> 
> pymc3_trace = pm.sample(1000, tune=1000, chains=4)
> pm.traceplot(pymc3_trace)
> 
> ```

---

<div class="post-metadata">

**Author:** ![Brian](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/brian/32/7088_2.png) [@Brian](https://discourse.pymc.io/u/Brian)\
**Post date:** [February 6, 2019, 1:46pm UTC](https://discourse.pymc.io/t/zero-inflated-poisson-log-lik/2664/3 "2019-02-06T13:46:47Z")

</div>

Thank you so much! I appreciate it!

One issue I also found was confusion on my side in regards to what ‘p’ represents. I was under the impression it was the probability of zero, but it seems as parameterized I need to pass in (1-p) to match what I did with the built in function (to match the rethinking statistics book) as this set-up in pymc3 appears to assume p is the probability of not zero?

---

<div class="post-metadata">

**Author:** ![katharina](https://avatars.discourse-cdn.com/v4/letter/k/7993a0/32.png) [@katharina](https://discourse.pymc.io/u/katharina)\
**Post date:** [August 2, 2019, 10:11am UTC](https://discourse.pymc.io/t/zero-inflated-poisson-log-lik/2664/4 "2019-08-02T10:11:47Z")

</div>

I stumbled over the p / psi parameters as well, maybe this can be helpful to someone (first post in a forum, please correct me if I am wrong).

In the [pxmc3 docs](https://docs.pymc.io/api/distributions/discrete.html#pymc3.distributions.discrete.ZeroInflatedPoisson) it is explained:  
“psi : float, Expected proportion of **_Poisson_** variates (0 \< psi \< 1)”  
In [Statistical Rethinking](https://speakerdeck.com/rmcelreath/l13-statistical-rethinking-winter-2019?slide=31) (and also the omega parameter in statsmodels. distributions.discrete.zipoisson), the p value refers to the probability of excess zeros generated by the **_Binomial_** proportion of the mixed model.

So if the probability p that excess zeros come from the binomial process is p=0.2, the probability that an observed zero came from the poisson process is 1-p=0.8.  
With theta (lambda) = 1, the probability of observing a zero overall is just the same in both formulations:  
PYMC3 formulation: (1-psi) + psi exp(-theta) = (1-0.8) + 0.8 exp(-1) = 0.4943035529371539  
Statistical Rethinking formulation: p + (1-p) exp(-lambda) = 0.2 + (1-0.2) exp(-1) = 0.4943035529371539