sample_prior_predictive not working with Negative Binomial likelihood

[AlexAndorra]

Description of your problem

As discussed here with @lucianopaz, the following model samples smoothly when using pm.sample and pm.sample_posterior_predictive. But it raises ValueError: operands could not be broadcast together with shapes (500,10886) (500,) when trying pm.sample_prior_predictive (10886 is the number of data points):

with pm.Model() as m_bike_poisson:
    a = pm.Normal("a", 0.0, 0.5)
    bT = pm.Normal("bT", 0.0, 0.2)

    lam = pm.math.exp(a + bT * bike_data["temp_std"])
    scale = pm.Exponential("scale", 2.0)

    bike_count = pm.NegativeBinomial(
        "bike_count", mu=lam, alpha=scale, observed=bike_data["count"]
    )

    prior_checks = pm.sample_prior_predictive(random_seed=RANDOM_SEED)
    trace_bike_poisson = pm.sample(1000, tune=2000, random_seed=RANDOM_SEED)
    post_samples = pm.sample_posterior_predictive(
        trace_bike_poisson, random_seed=RANDOM_SEED
    )

Adding shape=bike_data["temp_std"].to_numpy().shape to the Negative Binomial doesn't help.
The data simply come from Kaggle’s bike-sharing demand contest.
It looks like a shape issue when drawing random values from the Gamma distribution:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-3e62e449f5b3> in <module>
     16     )
     17 
---> 18     prior_checks = pm.sample_prior_predictive(random_seed=RANDOM_SEED)
     19     trace_bike_poisson = pm.sample(1000, tune=2000, random_seed=RANDOM_SEED)
     20     post_samples = pm.sample_posterior_predictive(

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1770     names = get_default_varnames(vars_, include_transformed=False)
   1771     # draw_values fails with auto-transformed variables. transform them later!
-> 1772     values = draw_values([model[name] for name in names], size=samples)
   1773 
   1774     data = {k: v for k, v in zip(names, values)}

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    632                                         point=point,
    633                                         givens=temp_givens,
--> 634                                         size=size)
    635                     givens[next_.name] = (next_, value)
    636                     drawn[(next_, size)] = value

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    831                         else:
    832                             dist_tmp.shape = val.shape
--> 833                 return dist_tmp.random(point=point, size=size)
    834             else:
    835                 return param.distribution.random(point=point, size=size)

~/opt/anaconda3/envs/stat-rethink-pymc3/lib/python3.7/site-packages/pymc3/distributions/discrete.py in random(self, point, size)
    674         """
    675         mu, alpha = draw_values([self.mu, self.alpha], point=point, size=size)
--> 676         g = generate_samples(stats.gamma.rvs, alpha, scale=mu / alpha,
    677                              dist_shape=self.shape,
    678                              size=size)

ValueError: operands could not be broadcast together with shapes (500,10886) (500,) 

Versions and main components

• PyMC3 Version: Master
• Theano Version: 1.0.4
• Python Version: 3.7.6
• Operating system: OS X
• How did you install PyMC3: pip

[fearghuskeeble]

Yes! I found this problem too once.

It's trying to do the division mu / alpha before it sorts the shapes out. The scipy function needs a wrapper like the _random method in, e.g. ZeroInflatedNegativeBinomial.

I have this, I will try to make a pull request for it.

[AlexAndorra]

Closing, as this is fixed now