Unable to draw prior samples from model

[mschmidt87]

If you have questions about a specific use case, or you are not sure whether this is a bug or not, please post it to our discourse channel: https://discourse.pymc.io

Description of your problem

I have problems with drawing prior samples from my model. The model I am using is the linear trend component of a time series model (Facebook’s prophet model). I am using the implementation in PyMC3 of the timeseers package authored by @MBrouns. I’d like to draw prior sample from the model, but I keep running into some Theano error that I can’t figure out.

Please provide a minimal, self-contained, and reproducible example.
Here is a minimal reproducer of the model (I basically stripped down the timeseers implementation to a minimal version so that it’s easier to understand here):

def linear_trend_model(X, y):
    t = X["t"].values
    n_changepoints = 1
    s = np.linspace(0, np.max(t), n_changepoints + 2)[1:-1]

    with pm.Model() as model:
        A = (t[:, None] > s) * 1.0

        delta = pm.Laplace("delta", 0, 0.05, shape=(n_changepoints,))
        k = pm.Normal("k", 0, 1.0)
        m = pm.Normal("m", 0, 5)
        gamma = -s * delta

        mu = (k + pm.math.sum(A * delta, axis=1)) * t + (m + pm.math.sum(A * gamma, axis=1))
        sigma = pm.HalfCauchy("sigma", 0.5)
        pm.Normal("obs", mu=mu, sd=sigma, observed=y)
    return model

To draw prior samples, I am creating some dummy data that I can feed into the model and then call pm.draw_sample_prior_predictive:

import pandas as pd
import numpy as np
import pymc3 as pm

t = np.linspace(0, 1, 24)
y = np.random.rand(len(t))
df = pd.DataFrame({'t': t, 'y': y})

with linear_trend_model(df, df.y):
    prior_trace = pm.sample_prior_predictive(10)

Please provide the full traceback.
This yields this error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    834                 try:
--> 835                     return dist_tmp.random(point=point, size=size)
    836                 except (ValueError, TypeError):

~//path/to/lib/python3.7/site-packages/pymc3/distributions/continuous.py in random(self, point, size)
    513         mu, tau, _ = draw_values([self.mu, self.tau, self.sigma],
--> 514                                  point=point, size=size)
    515         return generate_samples(stats.norm.rvs, loc=mu, scale=tau**-0.5,

~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    697                                             givens=givens.values(),
--> 698                                             size=size)
    699                         evaluated[param_idx] = drawn[(param, size)] = value

~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    875             func = _compile_theano_function(param, input_vars)
--> 876             output = func(*input_vals)
    877             return output

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
   2107 
-> 2108         return self._vectorize_call(func=func, args=vargs)
   2109 

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
   2181         if self.signature is not None:
-> 2182             res = self._vectorize_call_with_signature(func, args)
   2183         elif not args:

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call_with_signature(self, func, args)
   2222         for index in np.ndindex(*broadcast_shape):
-> 2223             results = func(*(arg[index] for arg in args))
   2224 

~//path/to/lib/python3.7/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    812                             arg, strict=s.strict,
--> 813                             allow_downcast=s.allow_downcast)
    814 

~//path/to/lib/python3.7/site-packages/theano/tensor/type.py in filter(self, data, strict, allow_downcast)
    194                 raise TypeError("Non-unit value on shape on a broadcastable"
--> 195                                 " dimension.", data.shape, self.broadcastable)
    196             i += 1

TypeError: ('Bad input argument to theano function with name "/Users/maximilianschmidt2//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py:732" at index 0 (0-based). ', 'Non-unit value on shape on a broadcastable dimension.', (10,), (True,))

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-11-4a0c0bbadcf5> in <module>
      1 with linear_trend_model(df, df.y):
----> 2     prior_trace = pm.sample_prior_predictive(10)

~//path/to/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1955     names = get_default_varnames(vars_, include_transformed=False)
   1956     # draw_values fails with auto-transformed variables. transform them later!
-> 1957     values = draw_values([model[name] for name in names], size=samples)
   1958 
   1959     data = {k: v for k, v in zip(names, values)}

~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    652                                         point=point,
    653                                         givens=temp_givens,
--> 654                                         size=size)
    655                     givens[next_.name] = (next_, value)
    656                     drawn[(next_, size)] = value

~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    842                     with _DrawValuesContextBlocker():
    843                         val = np.atleast_1d(dist_tmp.random(point=point,
--> 844                                                             size=None))
    845                     # Sometimes point may change the size of val but not the
    846                     # distribution's shape

~//path/to/lib/python3.7/site-packages/pymc3/distributions/continuous.py in random(self, point, size)
    512         """
    513         mu, tau, _ = draw_values([self.mu, self.tau, self.sigma],
--> 514                                  point=point, size=size)
    515         return generate_samples(stats.norm.rvs, loc=mu, scale=tau**-0.5,
    516                                 dist_shape=self.shape,

~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    696                                             point=point,
    697                                             givens=givens.values(),
--> 698                                             size=size)
    699                         evaluated[param_idx] = drawn[(param, size)] = value
    700                         givens[param.name] = (param, value)

~//path/to/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    874                 input_vals = []
    875             func = _compile_theano_function(param, input_vars)
--> 876             output = func(*input_vals)
    877             return output
    878     raise ValueError('Unexpected type in draw_value: %s' % type(param))

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
   2106             vargs.extend([kwargs[_n] for _n in names])
   2107 
-> 2108         return self._vectorize_call(func=func, args=vargs)
   2109 
   2110     def _get_ufunc_and_otypes(self, func, args):

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
   2180         """Vectorized call to `func` over positional `args`."""
   2181         if self.signature is not None:
-> 2182             res = self._vectorize_call_with_signature(func, args)
   2183         elif not args:
   2184             res = func()

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call_with_signature(self, func, args)
   2210 
   2211         broadcast_shape, dim_sizes = _parse_input_dimensions(
-> 2212             args, input_core_dims)
   2213         input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
   2214                                          input_core_dims)

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in _parse_input_dimensions(args, input_core_dims)
   1871     dim_sizes = {}
   1872     for arg, core_dims in zip(args, input_core_dims):
-> 1873         _update_dim_sizes(dim_sizes, arg, core_dims)
   1874         ndim = arg.ndim - len(core_dims)
   1875         dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])

~//path/to/lib/python3.7/site-packages/numpy/lib/function_base.py in _update_dim_sizes(dim_sizes, arg, core_dims)
   1837             '%d-dimensional argument does not have enough '
   1838             'dimensions for all core dimensions %r'
-> 1839             % (arg.ndim, core_dims))
   1840 
   1841     core_shape = arg.shape[-num_core_dims:]

ValueError: 0-dimensional argument does not have enough dimensions for all core dimensions ('i_1_0',)

Please provide any additional information below.

Versions and main components

• PyMC3 Version: 3.9.3
• Theano Version: 1.0.5
• Python Version: 3.7.7
• Operating system: Mac OS
• How did you install PyMC3: pip

[junpenglao]

Sounds like this is related to the bug of how we are treating (1, ) tensor. Maybe #4214 would fix it? @michaelosthege

[MarcoGorelli]

works for n_changepoints > 1

[michaelosthege]

Yes, #4214 should fix it.

It's possible to do a workaround, but very nasty.. I faced the same problem today and decided to work with the branch from #4214. It will be one of the first things that we merge after the 3.10 release is out.

[mschmidt87]

Thank you for your replies. I could confirm that:

• my minimal reproducer works for n_changepoints > 1 (as found by @MarcoGorelli )
• with Fixing scalar shape handling: treat (1,) as vectors #4214 , it works for n_changepoints = 1.
• Also, using the full implementation from the timeseers package, I can draw prior samples from the model using Fixing scalar shape handling: treat (1,) as vectors #4214 .

[michaelosthege]

I'm closing this because it was already fixed on master. The fix will be released with 3.11.0.