Unable to draw prior samples from model

Questions
1

Hi,

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 timeseerspackage 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.

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)

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',)
2

Yes, I can confirm this problem. Please file a bug report on our issue tracker

3

Thank you for your fast reply. I filed an issue here: https://github.com/pymc-devs/pymc3/issues/4253

4

I’ve only looked at the example briefly, but I think I know what the problem is.
In order to keep the implementation in timeseers simple I always assume there is an additional dimension on the parameters for estimating unpooled effects. Note the shape=(n_groups, ...) on for example the definition of delta.

delta = pm.Laplace(
    self._param_name("delta"), 0, self.changepoints_prior_scale, shape=(n_groups, self.n_changepoints)
)
k = pm.Normal(self._param_name("k"), 0, self.growth_prior_scale, shape=n_groups)

Because of this, I don’t use regular dot products in the definition of mu, but I use pm.math.sum(A * delta). If you don’t want this extra dimension on your parameters, you should probably change your definition of mu accordingly. e.g.:

g = (k + dot(A, delta)) * t + (m + dot(A, gamma))
5

Thank you for your reply. In my minimal example, I removed the groups, however, ultimately I’d like to keep that dimension because I do want to use partial pooling, so just removing that extra dimension is not really what I’m looking for. And I am getting the same error when using the timeseers implementation.

I also don’t quite understand your comment: In the implementation in timeseers (https://github.com/MBrouns/timeseers/blob/50275a4f50eae9b7af0083b8db91d57b63684692/src/timeseers/linear_trend.py#L48), for instance, A.shape = (len(t), n_changepoints), right? And delta.shape = (n_groups, n_changepoints), so that delta[group].shape = (len(t), n_changepoints). So if you want to sum out the n_changepoints dimension, using the sum(A, delta[group], axis=1) approach is the only way I can think of. But perhaps I misunderstood what you tried to say.

6

You’re right, I should’ve looked more carefully! I didn’t realize that sample_prior_predictive didn’t work while pm.sample did work. I’m not super familiar with that piece of the PyMC codebase, but couldn’t you, in that case, define the model without any observed= variable and use pm.sample as normal?

7

As it was pointed out in the issue on github (https://github.com/pymc-devs/pymc3/issues/4253), this issue is fixed completely (as far as I can see) by this PR: https://github.com/pymc-devs/pymc3/pull/4214 .

That is, both this minimal reproducer as well as the full implementation in timeseers work using the branch from the PR.

1 Like