Zero Inflated Poisson Log-Lik


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

#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?