TypeError on sample_ppc with binomial random variable

[kzielnicki]

Minimal example, simple "coin flip" model:

import pymc3 as pm
import numpy as np
import scipy.stats as ss

data = ss.binom(1, 0.8).rvs(20)

with pm.Model() as model:
    # prior for p
    p = pm.Uniform('p', lower=0, upper=1)
    
    # model for x given p
    x = pm.Binomial('x', n=np.ones(len(data)), p=p, observed=data)

with model:
    trace = pm.sample(draws=1000, init='advi', n_init=20000)

ppc = pm.sample_ppc(trace, model=model, samples=500, size=100)

Gives:

TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'

Full trace:

--------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-46-d7601fc9ec00> in <module>()
----> 1 ppc = pm.sample_ppc(trace, model=model)

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/pymc3/sampling.pyc in sample_ppc(trace, samples, model, vars, size, random_seed, progressbar)
    389         for var in vars:
    390             ppc[var.name].append(var.distribution.random(point=param,
--> 391                                                          size=size))
    392 
    393     return {k: np.asarray(v) for k, v in ppc.items()}

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/pymc3/distributions/discrete.pyc in random(self, point, size, repeat)
     49         return generate_samples(stats.binom.rvs, n=n, p=p,
     50                                 dist_shape=self.shape,
---> 51                                 size=size)
     52 
     53     def logp(self, value):

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in generate_samples(generator, *args, **kwargs)
    349                                         [int(x) for x in broadcast_shape],
    350                                         [int(x) for x in prefix_shape],
--> 351                                         *args, **kwargs)
    352             if broadcast_shape == (1,):
    353                 samples = np.reshape(samples, prefix_shape)

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in replicate_samples(generator, size, repeats, *args, **kwargs)
    269     n = int(np.prod(repeats))
    270     if n == 1:
--> 271         samples = generator(size=size, *args, **kwargs)
    272     else:
    273         samples = np.array([generator(size=size, *args, **kwargs)

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.pyc in rvs(self, *args, **kwargs)
   2808         """
   2809         kwargs['discrete'] = True
-> 2810         return super(rv_discrete, self).rvs(*args, **kwargs)
   2811 
   2812     def pmf(self, k, *args, **kwds):

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.pyc in rvs(self, *args, **kwds)
    952         # by _rvs().
    953         self._size = size
--> 954         vals = self._rvs(*args)
    955 
    956         vals = vals * scale + loc

/Users/kevinzielnicki/python/virtualenv/stitch/lib/python2.7/site-packages/scipy/stats/_discrete_distns.pyc in _rvs(self, n, p)
     38     """
     39     def _rvs(self, n, p):
---> 40         return self._random_state.binomial(n, p, self._size)
     41 
     42     def _argcheck(self, n, p):

mtrand.pyx in mtrand.RandomState.binomial (numpy/random/mtrand/mtrand.c:31765)()

TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'

[aseyboldt]

Yes, I think this is a bug. You can work around it by changing the type of n:

    # model for x given mu
    x = pm.Binomial('x', n=np.ones(len(data), dtype='int'),
                    p=p, observed=data, shape=20)

I'll have a closer look next week if nobody beats me to it.

[aseyboldt]

@junpenglao I don't think this should be a blocker for 3.1. The current behaviour is not really that unreasonable, the error message could be better though.

This raises the question how we want to deal with real values for n in binomial. The term used for logp spits out values for real n, as it uses the gamma function internally. Strictly speaking, this doesn't result in a normalized distribution for non-integer n. So we could just check if n is of integer type and raise an error if it is not. But using non-integer values for n can be useful, as this allows us to use NUTS to estimate it and if n is sufficiently large and p isn't close to 1, the error should be quite small.

I guess we could add a new keyword allow_real_n, and if that is not set to True, we throw an error for bad n?

[junpenglao]

Hmm, I would suggest converting n to integer internally, which is this PR is currently working on.

And yeah so this should not be a blocker for 3.1.

[aseyboldt]

I don't think we should go around just casting values. If the type supplied by a user is not what we need, we should throw an error.
Also, if we cast n to an integer, we can't use nuts anymore.

[junpenglao]

Also, if we cast n to an integer, we can't use nuts anymore.

I am not sure that would be the case, if you just use n=10 nuts can still runs

[aseyboldt]

But the gradients are zero. The derivative of floor is zero almost everywhere. That completely breaks nuts.

[junpenglao]

Oh you mean if n being an RV in the graph - yeah in that case we dont want to cast it into int, hmm I think I get what you mean now...

[junpenglao]

I wonder what is the solution in Stan in these cases.

[aseyboldt]

Stan simply doesn't allow n to be anything but an integer.