Custom discrete distribution for zero truncated Poisson

Questions
1

I am trying to create a custom distribution for zero truncated Poisson by specifying the log probability function and then using the pm.DensityDist method as below. However I get these errors:

TypeError: don't know how to convert scalar number to int

k = np.array([1,2, 48,1])
def log_fact(a):
    return np.array([ np.sum([np.log(i) for i in range(1, x)]) for x in np.nditer(a, ['refs_ok']) ]).reshape(a.shape)

print(log_fact(k))

lam = 1
def ztp_logp(k):
    return ( k * np.log(lam) - np.log(np.exp(lam) - 1) - log_fact(k) ).sum()

ztp_logp(k)


with pm.Model() as ztp_model:
    zp_poisson = pm.DensityDist('zp_poisson', ztp_logp, observed=k)

TypeErrorTraceback (most recent call last)
<ipython-input-66-43b7242eea0e> in <module>()
      1 with pm.Model() as ztp_model:
----> 2     zp_poisson = pm.DensityDist('zp_poisson', ztp_logp, observed=k)

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
     35             total_size = kwargs.pop('total_size', None)
     36             dist = cls.dist(*args, **kwargs)
---> 37             return model.Var(name, dist, data, total_size)
     38         else:
     39             raise TypeError("Name needs to be a string but got: {}".format(name))

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/model.pyc in Var(self, name, dist, data, total_size)
    780                 var = ObservedRV(name=name, data=data,
    781                                  distribution=dist,
--> 782                                  total_size=total_size, model=self)
    783             self.observed_RVs.append(var)
    784             if var.missing_values:

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/model.pyc in __init__(self, type, owner, index, name, data, distribution, total_size, model)
   1228 
   1229             self.missing_values = data.missing_values
-> 1230             self.logp_elemwiset = distribution.logp(data)
   1231             # The logp might need scaling in minibatches.
   1232             # This is done in `Factor`.

<ipython-input-64-5d9727a596ae> in ztp_logp(k)
      4 from scipy.special import factorial
      5 def ztp_logp(k):
----> 6     return ( k * np.log(lam) - np.log(np.exp(lam) - 1) - log_fact(k) ).sum()
      7 
      8 ztp_logp(k)

<ipython-input-63-03136e096fc4> in log_fact(a)
      1 k = np.array([1,2, 48,1])
      2 def log_fact(a):
----> 3     return np.array([ np.sum([np.log(i) for i in range(1, x)]) for x in np.nditer(a, ['refs_ok']) ]).reshape(a.shape)
      4 
      5 print(log_fact(k))

TypeError: don't know how to convert scalar number to int
1 Like
2

Your log_fact implementation is likely not very efficient. Try instead use gamma function instead of factorial. You should also try to rewriting it into a vectorised numpy version, and then replace the numpy operation with theano operation.

3

This should work:

import theano.tensor as tt
def ztp_logp(k):
    return tt.sum(k * tt.log(lam) - tt.log(tt.exp(lam) - 1) - tt.gammaln(k))
4

Hi Junpenglao,
Thanks for the quick response. Using theano tensor and the gamma log appeared to solve that specific issue with conversion from scalar to int.
However, now it complains at the posterior sampling step. I am assuming for this Zero Truncated Poisson distribution that I am intending to create, there is something missing?

k = np.array([1,2, 48,1])
import theano.tensor as tt
def ztp_logp(k):
    return tt.sum(k * tt.log(lam) - tt.log(tt.exp(lam) - 1) - tt.gammaln(k))

ztp_logp(k)

Sum{acc_dtype=float64}.0

with pm.Model() as ztp_model:
    zp_poisson = pm.DensityDist('zp_poisson', ztp_logp, observed=k)
    zstrace = pm.sample(5000)
    pm.plot_posterior(zstrace)

ValueErrorTraceback (most recent call last)
<ipython-input-96-e13a83ea7cad> in <module>()
      1 with pm.Model() as ztp_model:
      2     zp_poisson = pm.DensityDist('zp_poisson', ztp_logp, observed=k)
----> 3     zstrace = pm.sample(5000)
      4     pm.plot_posterior(zstrace)

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/sampling.pyc in sample(draws, step, init, n_init, start, trace, chain, njobs, tune, nuts_kwargs, step_kwargs, progressbar, model, random_seed, live_plot, discard_tuned_samples, live_plot_kwargs, **kwargs)
    242 
    243     if model.ndim == 0:
--> 244         raise ValueError('The model does not contain any free variables.')
    245 
    246     if step is None and init is not None and pm.model.all_continuous(model.vars):

ValueError: The model does not contain any free variables.
5

As the error message suggest, you need to include free variables for the sampler to do inference on. For example, a valid model would look something like this:

with pm.Model():
    lam = pm.HalfNormal('lam', sd=10.)
    def ztp_logp(k):
        return tt.sum(k * tt.log(lam) - tt.log(tt.exp(lam) - 1) - tt.gammaln(k))
    like = pm.DensityDist('ztp', ztp_logp, observed=k)
    trace = pm.sample()

where lam is a free parameter. It was fixed to one in your case above which cause the error.
For more information about how to use Densitydist, you can do a search in the doc or discourse.

6

Excellent, that is very helpful. I looked for examples in the docs, but could not find one explained all the way.
I could add this as an example to the docs, showing the creation for a discrete random variable (zero truncated poisson). Now the next step that I am struggling with is specifying a random function, in order to perform posterior predictive checks (and drawing samples).
This is what I have so far:

import theano.tensor as tt
class ZeroTruncatedPoisson(pm.Discrete):
    '''
    Similar to a Poisson distribution, but the support is over positive integes 
    and so, excludes zero. An example, number of items in your grocery shopping cart
    at checkout
    '''
    def __init__(self, lam):
        self._lam = tt.as_int_none_variable(lam)

    def logp(self, value):
        '''
        compute total_log_probability
        '''
        return tt.sum(value * tt.log(self._lam) - tt.log(tt.exp(self._lam) -1) - tt.gammaln(value))
    
    def random(self, point=None, size=None, repeat=None):
        '''
        Generates a random sample from Zero Truncated Poisson
        '''
        k = 1
        t = tt.exp(-self._lam)/(1 - tt.exp(-self._lam))/(1. * self._lam)
        s = t
        u = np.random.uniform()
        while s < u:
            k = k + 1
            t = t * self._lam / k
            s = s + t
            u = np.random.uniform()
        return k
    

with pm.Model() as zmodel:
    zp_poisson = ZeroTruncatedPoisson('zp', lam=1, observed=k)

But this fails with:


AttributeErrorTraceback (most recent call last)
<ipython-input-118-2be01f1fdb32> in <module>()
      1 with pm.Model() as zmodel:
----> 2     zp_poisson = ZeroTruncatedPoisson('zp', lam=1, observed=k)

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
     35             total_size = kwargs.pop('total_size', None)
     36             dist = cls.dist(*args, **kwargs)
---> 37             return model.Var(name, dist, data, total_size)
     38         else:
     39             raise TypeError("Name needs to be a string but got: {}".format(name))

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/model.pyc in Var(self, name, dist, data, total_size)
    780                 var = ObservedRV(name=name, data=data,
    781                                  distribution=dist,
--> 782                                  total_size=total_size, model=self)
    783             self.observed_RVs.append(var)
    784             if var.missing_values:

/Users/vmullachery/anaconda/envs/dl2.7/lib/python2.7/site-packages/pymc3/model.pyc in __init__(self, type, owner, index, name, data, distribution, total_size, model)
   1218         if type is None:
   1219             data = pandas_to_array(data)
-> 1220             type = TensorType(distribution.dtype, data.shape)
   1221 
   1222         self.observations = data

AttributeError: 'ZeroTruncatedPoisson' object has no attribute 'dtype'
7

Dear junpenglao,
What causes the previous error?

AttributeError: 'ZeroTruncatedPoisson' object has no attribute 'dtype'

I am trying to create a class that has logp and random methods for ZeroTruncatedPoisson

8

You can add a **kwargs in the init and let the default take care of it:

    def __init__(self, lam, **kwargs):
        self._lam = tt.as_int_none_variable(lam)
        super(ZeroTruncatedPoisson, self).__init__(**kwargs)
9

@junpenglao, thank you.
Finally I was able to get a ZeroTruncatedPoisson working for my requirements (not generalized to multiple dimensions):

import theano.tensor as tt
from theano import function
from theano import printing
import pymc3 as pm

class ZeroTruncatedPoisson(pm.Discrete):
    '''
    Similar to a Poisson distribution, but the support is over positive integes 
    and so, excludes zero. An example, number of items in your grocery shopping cart
    at checkout
    '''
    def __init__(self, mu, *args, **kwargs):
        super(ZeroTruncatedPoisson, self).__init__(*args, **kwargs)
        self.mu = tt.as_tensor_variable(mu)
    
    def logp(self, value):
        '''
        compute total_log_probability
        '''
        return tt.sum(value * tt.log(self.mu) - tt.log(tt.exp(self.mu) -1) - tt.gammaln(value))
    
    def random(self, point=None, size=None):
        '''
        Generates a random sample from Zero Truncated Poisson
        '''
        k = 1
        mu = self.mu.tag.test_value
        t = tt.exp(-self.mu)/(1 - tt.exp(-self.mu))/(1. * self.mu)
        tn = t.tag.test_value
        s = tn
        u = np.random.uniform()
        while s < u:
            k = k + 1
            tn = tn * mu / k
            s = s + tn
            u = np.random.uniform()
        return k

And then infer the posterior and generate ppc samples:

obs = np.array([1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,48,1])
with pm.Model() as zmodel:
    mu = pm.HalfNormal('mu', sd=10.)
    zp_poisson = ZeroTruncatedPoisson('zp', mu=mu, observed=obs)
    ztrace = pm.sample(1000)
    pm.plot_posterior(ztrace)

with zmodel:
    ppc = pm.sample_ppc(ztrace, 5000)
    sns.distplot(ppc['zp'], label='Simulated')
    sns.distplot(obs, label='Real')
    plt.legend()
1 Like
10

Why does NUTS work without the gradient function grad() in class ZeroTruncatedPoisson? It was required to write grad() in the past releases of pymc3. What has changed?

11

If you write the logp in theano (see the line return tt.sum(value * tt.log(self.mu) - tt.log(tt.exp(self.mu) -1) - tt.gammaln(value))), theano will take the gradient automatically so you dont need to implement custom grad() method.

1 Like