Behavior of GaussianRandomWalk with tensor parameters

[AlphaScorpii-dev]

I was trying to look after some trend changes in random walks, using a model similar to the one below.
When using a Normal RV for data, things work fine. The break location and the different mu/sigma are estimated properly.
Now, when using the random walk and observing the cumsum values, everything breaks.
Somehow, it seems that GaussianRandowWalk returns matrices instead of vectors during a run, when passing tensors as mu and sigma.
I assume I am missing something obvious. I tried using pandas.Series as observed data as well, to no avail.

Please provide a minimal, self-contained, and reproducible example.

import numpy, pymc

with pymc.Model() as model:
    ...:     n=1000
    ...:     d=numpy.random.randn(n)
    ...:     m1=pymc.Normal("m1", mu=0, sigma=1)
    ...:     m2=pymc.Normal("m2", mu=0, sigma=1)
    ...:     s1=pymc.Exponential("s1", lam=1)
    ...:     s2=pymc.Exponential("s2", lam=1)
    ...:     switchpoint=pymc.DiscreteUniform("switchpoint", lower=n/10, upper=n*9/10)
    ...:     m=pymc.math.switch(switchpoint>numpy.arange(n), m1, m2)
    ...:     s=pymc.math.switch(switchpoint>numpy.arange(n), s1, s2)
    ...:     #works nicely with gaussian RV
    ...:     #data=pymc.Normal("data", mu=m, sigma=s, observed=d)
    ...:     data=pymc.GaussianRandomWalk("data", mu=m, sigma=s, observed=numpy.cumsum(d))
    ...:     trace=pymc.sample(10000)

Please provide the full traceback.

Complete error traceback

C:\Users\me\AppData\Local\miniforge3\envs\science\lib\site-packages\pymc\distributions\timeseries.py:275: UserWarning: Initial distribution not specified, defaulting to `Normal.dist(0, 100)`.You can specify an init_dist manually to suppress this warning.
  warnings.warn(
---------------------------------------------------------------------------
ShapeError                                Traceback (most recent call last)
Input In [81], in <cell line: 1>()
      7 m=pymc.math.switch(switchpoint>numpy.arange(n), m1, m2)
      8 s=pymc.math.switch(switchpoint>numpy.arange(n), s1, s2)
----> 9 data=pymc.GaussianRandomWalk("data", mu=m, sigma=s, observed=d)
     10 trace=pymc.sample(10000)

File ~\AppData\Local\miniforge3\envs\science\lib\site-packages\pymc\distributions\timeseries.py:249, in GaussianRandomWalk.__new__(cls, steps, *args, **kwargs)
    241 def __new__(cls, *args, steps=None, **kwargs):
    242     steps = get_steps(
    243         steps=steps,
    244         shape=None,  # Shape will be checked in `cls.dist`
   (...)
    247         step_shape_offset=1,
    248     )
--> 249     return super().__new__(cls, *args, steps=steps, **kwargs)

File ~\AppData\Local\miniforge3\envs\science\lib\site-packages\pymc\distributions\distribution.py:271, in Distribution.__new__(cls, name, rng, dims, initval, observed, total_size, transform, *args, **kwargs)
    267 if resize_shape:
    268     # A batch size was specified through `dims`, or implied by `observed`.
    269     rv_out = change_rv_size(rv=rv_out, new_size=resize_shape, expand=True)
--> 271 rv_out = model.register_rv(
    272     rv_out,
    273     name,
    274     observed,
    275     total_size,
    276     dims=dims,
    277     transform=transform,
    278     initval=initval,
    279 )
    281 # add in pretty-printing support
    282 rv_out.str_repr = types.MethodType(str_for_dist, rv_out)

File ~\AppData\Local\miniforge3\envs\science\lib\site-packages\pymc\model.py:1379, in Model.register_rv(self, rv_var, name, data, total_size, dims, transform, initval)
   1372         raise TypeError(
   1373             "Variables that depend on other nodes cannot be used for observed data."
   1374             f"The data variable was: {data}"
   1375         )
   1377     # `rv_var` is potentially changed by `make_obs_var`,
   1378     # for example into a new graph for imputation of missing data.
-> 1379     rv_var = self.make_obs_var(rv_var, data, dims, transform)
   1381 return rv_var

File ~\AppData\Local\miniforge3\envs\science\lib\site-packages\pymc\model.py:1405, in Model.make_obs_var(self, rv_var, data, dims, transform)
   1402 data = convert_observed_data(data).astype(rv_var.dtype)
   1404 if data.ndim != rv_var.ndim:
-> 1405     raise ShapeError(
   1406         "Dimensionality of data and RV don't match.", actual=data.ndim, expected=rv_var.ndim
   1407     )
   1409 if aesara.config.compute_test_value != "off":
   1410     test_value = getattr(rv_var.tag, "test_value", None)

ShapeError: Dimensionality of data and RV don't match. (actual 1 != expected 2)

Please provide any additional information below.

Versions and main components

• PyMC/PyMC3 Version: 4.1.2
• Aesara/Theano Version: 2.7.5
• Python Version: 3.10.5
• Operating system: Win10 21H2
• How did you install PyMC/PyMC3: (conda/pip) conda & conda-forge repository

[ricardoV94]

The GRW expects a single mu/sigma per independent series.

[ricardoV94]

If you want to model a RW with changing mu you will have to do it manually. If you difference your data, you can use a normal likelihood.

[AlphaScorpii-dev]

Thanks for your feedback @ricardoV94
If you allow me a question then, what is the fundamental difference between Normal and GaussianRandomWalk, that would enable one to use series, and the other to be limited to scalars only ?
Looking at the source for GaussianRandomWalkRV, the class has attributes that seem to limit the dimension

ndim_supp = 1
ndims_params = [0, 0, 0, 0]

where other RV like the ones provided in aesara would rather be like

ndim_supp = 0
ndims_params = [0, 0]

Is this how I can spot the difference ?

[ricardoV94]

Yes, that is where the difference comes from. It was a simplifying choice to make the distribution be a vector one (ndim_supp==1) with shared scalar mu and sigma (ndims_params=[0, 0, ...])

Your interpretation of the GRW is equally valid, and perhaps more general. But it would make it more hazardous to write down batched GRW, as the user would need to know/remember to broadcast mu and sigma along the last axis if they didn't mean it to vary within the same series (if ndim_supp==0), which I think is the most common case. Or alternatively (if ndims_params=[1, 1, ...]), they would always need to pass a vector for mu and sigma.

This guide explains the distribution dimensionality logic: https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/dimensionality.html

Anyway, you are free to adapt the source code and implement the version that matches your interpretation. If other users feel the need for it, we can even consider adding it to the codebase or replace the current one altogether.

[ricardoV94]

Follow up discussion in #6090