Sample_prior_predictive doesn't work with Negative Binomial likelihood

Hi all,
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
    )

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, but I wanted to be sure before raising a GitHub issue – I think @lucianopaz is the go-to person here?
Thanks in advance :vulcan_salute:

Does it also fail if you explicitly pass the shape to the negative binomial? Something like shape=bike_data["temp_std"].to_numpy().shape?

Yes, same error I’m afraid :confused:
Do you want me to file a GH issue?

I think that would be the best thing to do. Although, with the covid quarantine and trying to deal with working from home, and all the family related stuff at the same time, I won’t be able to look into pymc development for some time. Hopefully someone else will be able to pitch in.

Just did it, thanks Luciano.
Will try and take a look, but same covid situation and working on other PRs.
Hope you and your family are safe, and good luck for the quarantine!