# Custom discrete distribution for zero truncated Poisson

**URL:** <https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622>\
**Category:** Questions\
**Created:** [December 2, 2017, 1:54pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622 "2017-12-02T13:54:26Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![mullachv](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mullachv/32/400_2.png) [@mullachv](https://discourse.pymc.io/u/mullachv)\
**Post date:** [December 2, 2017, 1:54pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/1 "2017-12-02T13:54:26Z")

</div>

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:

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

```

```auto
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)

```

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

```

```auto
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

```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 2, 2017, 5:56pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/2 "2017-12-02T17:56:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 2, 2017, 7:15pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/3 "2017-12-02T19:15:54Z")

</div>

This should work:

```python
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))

```

---

<div class="post-metadata">

**Author:** ![mullachv](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mullachv/32/400_2.png) [@mullachv](https://discourse.pymc.io/u/mullachv)\
**Post date:** [December 3, 2017, 3:24am UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/4 "2017-12-03T03:24:04Z")

</div>

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?

```auto
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)

```

```auto
Sum{acc_dtype=float64}.0

```

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

```

```auto
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.

```

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 3, 2017, 7:28am UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/5 "2017-12-03T07:28:38Z")

</div>

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:

```python
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.

---

<div class="post-metadata">

**Author:** ![mullachv](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mullachv/32/400_2.png) [@mullachv](https://discourse.pymc.io/u/mullachv)\
**Post date:** [December 3, 2017, 7:46am UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/6 "2017-12-03T07:46:46Z")

</div>

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:

```auto
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
    

```

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

```

But this fails with:

```auto

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'

```

---

<div class="post-metadata">

**Author:** ![mullachv](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mullachv/32/400_2.png) [@mullachv](https://discourse.pymc.io/u/mullachv)\
**Post date:** [December 3, 2017, 7:40pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/7 "2017-12-03T19:40:59Z")

</div>

Dear junpenglao,  
What causes the previous error?

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

```

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

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 5, 2017, 9:41am UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/8 "2017-12-05T09:41:32Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![mullachv](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/mullachv/32/400_2.png) [@mullachv](https://discourse.pymc.io/u/mullachv)\
**Post date:** [December 6, 2017, 12:01am UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/9 "2017-12-06T00:01:25Z")

</div>

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

```auto
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:

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![Vlad](https://avatars.discourse-cdn.com/v4/letter/v/3bc359/32.png) [@Vlad](https://discourse.pymc.io/u/Vlad)\
**Post date:** [December 19, 2017, 3:23pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/10 "2017-12-19T15:23:31Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 19, 2017, 3:43pm UTC](https://discourse.pymc.io/t/custom-discrete-distribution-for-zero-truncated-poisson/622/11 "2017-12-19T15:43:35Z")

</div>

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.