Zero Inflated Poisson Log-Lik

Questions
1

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)
2

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 
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)
1 Like
3

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?

4

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 it is explained:
“psi : float, Expected proportion of Poisson variates (0 < psi < 1)”
In Statistical Rethinking (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

4 Likes