Newly introduced bug in fast_sample_posterior_predictive()

[bridgeland]

Description of your problem

There is a new bug in fast_sample_posterior_predictive() in pymc version 3.11.2, a bug that was not present as recently as version 3.9.3. The bug is not present in sample_posterior_predictive().

Here’s a minimal example that exhibits the bug. Note that the example uses @lucianopaz 's model factory approach, with two somewhat different models, one sampled and one sampled posterior predictive:

import pymc3 as pm 

b_observed = 1
with pm.Model() as m: 
    p = pm.Beta('p', 2, 2)
    b = pm.Binomial('b', n=1, p=p, observed=b_observed)
    trace = pm.sample(
        draws=600, tune=1000, chains=2, return_inferencedata=True)

with pm.Model() as m_forward:
    p = pm.Beta('p', 2, 2) 
    b2 = pm.Binomial('b2', n=1, p=p)
    b3 = pm.Binomial('b3', n=1, p=p*b2)

with m_forward:
    trace_forward = pm.fast_sample_posterior_predictive(
        trace, var_names=['p', 'b2', 'b3'])

This call to fast_sample_posterior_predictive() results in an AssertionError:

Complete error traceback

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-4-cbc865a15399> in <module>
      1 with m_forward:
----> 2     trace_forward = pm.fast_sample_posterior_predictive(trace, var_names=['p', 'b2', 'b3'])

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in fast_sample_posterior_predictive(trace, samples, model, var_names, keep_size, random_seed)
    289             strace = _trace if s == len_trace else _trace[slice(0, s)]
    290             try:
--> 291                 values = posterior_predictive_draw_values(cast(List[Any], vars), strace, s)
    292                 new_trace: dict[str, np.ndarray] = {k.name: v for (k, v) in zip(vars, values)}
    293                 ppc_trace.extend_trace(new_trace)

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in posterior_predictive_draw_values(vars, trace, samples)
    305 ) -> list[np.ndarray]:
    306     with _PosteriorPredictiveSampler(vars, trace, samples, None) as sampler:
--> 307         return sampler.draw_values()
    308 
    309 

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in draw_values(self)
    408                         # This may fail for autotransformed RVs, which don't
    409                         # have the random method
--> 410                         value = self.draw_value(next_, trace=trace, givens=temp_givens)
    411                         assert isinstance(value, np.ndarray)
    412                         givens[next_.name] = (next_, value)

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in draw_value(self, param, trace, givens)
    581                 assert isinstance(model, Model)
    582                 shape: tuple[int, ...] = tuple(_param_shape(param, model))
--> 583                 return random_sample(param.random, param, point=trace, size=samples, shape=shape)
    584             elif (
    585                 hasattr(param, "distribution")

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in random_sample(meth, param, point, size, shape)
    550             shape: tuple[int, ...],
    551         ) -> np.ndarray:
--> 552             val = meth(point=point, size=size)
    553             try:
    554                 assert val.shape == (size,) + shape, (

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/model.py in __call__(self, *args, **kwargs)
    103 
    104     def __call__(self, *args, **kwargs):
--> 105         return getattr(self.obj, self.method_name)(*args, **kwargs)
    106 
    107 

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/discrete.py in random(self, point, size)
    122         array
    123         """
--> 124         n, p = draw_values([self.n, self.p], point=point, size=size)
    125         return generate_samples(stats.binom.rvs, n=n, p=p, dist_shape=self.shape, size=size)
    126 

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    702     if ppc_sampler is not None:
    703         # this is being done inside new, vectorized sample_posterior_predictive
--> 704         return ppc_sampler(params, trace=point, samples=size)
    705 
    706     if point is None:

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in posterior_predictive_draw_values(vars, trace, samples)
    305 ) -> list[np.ndarray]:
    306     with _PosteriorPredictiveSampler(vars, trace, samples, None) as sampler:
--> 307         return sampler.draw_values()
    308 
    309 

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in draw_values(self)
    408                         # This may fail for autotransformed RVs, which don't
    409                         # have the random method
--> 410                         value = self.draw_value(next_, trace=trace, givens=temp_givens)
    411                         assert isinstance(value, np.ndarray)
    412                         givens[next_.name] = (next_, value)

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in draw_value(self, param, trace, givens)
    581                 assert isinstance(model, Model)
    582                 shape: tuple[int, ...] = tuple(_param_shape(param, model))
--> 583                 return random_sample(param.random, param, point=trace, size=samples, shape=shape)
    584             elif (
    585                 hasattr(param, "distribution")

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in random_sample(meth, param, point, size, shape)
    564                     pass
    565                 else:
--> 566                     raise e
    567 
    568             return val

~/opt/miniconda3/envs/brussels2_j/lib/python3.9/site-packages/pymc3/distributions/posterior_predictive.py in random_sample(meth, param, point, size, shape)
    552             val = meth(point=point, size=size)
    553             try:
--> 554                 assert val.shape == (size,) + shape, (
    555                     "Sampling from random of %s yields wrong shape" % param
    556                 )

AssertionError: Sampling from random of b2 ~ Binomial yields wrong shape

The AssertionError is raised in this code, within pymc3/distributions/posterior_predictive.py:

    def random_sample(
        meth: Callable[..., np.ndarray],
        param,
        point: _TraceDict,
        size: int,
        shape: tuple[int, ...],
    ) -> np.ndarray:
        val = meth(point=point, size=size)
        try:
            assert val.shape == (size,) + shape, (
                "Sampling from random of %s yields wrong shape" % param
            )

The debugger reveals values for the key variables:

ipdb>  val.shape
(1200,)
ipdb>  size
(1200,)
ipdb>  shape
()
ipdb>  size + shape
(1200,)

I’m not sure whether size should be 1200 instead of (1200,), or whether
the assertion should be assert val.shape == size + shape instead of assert val.shape == (size,) + shape. But something is wrong.

Note that this error was not raised in v3.9.3, and is not raised in v3.11.2 when sample_posterior_predictive() is called instead.

Versions and main components

• PyMC3 Version: 3.11.2
• Aesara/Theano Version: 1.1.2
• Python Version: 3.9.4
• Operating system: macos 10.15.7
• How did you install PyMC3: conda, I think (p=0.8)

[michaelosthege]

Closed by #4927