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

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?