Bug in the Multivariate Normal new random method?

Questions
1

Hi,

I am using the 3.10.0 pymc3 version and I know that the random method for the MvNormal was modified. I was using it and I faced a somewhat strange behavior. If I don’t explicitly define the shape argument I get a shape error, despite the fact that both mu and cov have a defined shape that could be inferred to do the sampling. From what I can see this happens because dist_shape is set to be shapeless when the shape argument is not supplied.

Example,

sigma = 2.0
cov = pm.gp.cov.WhiteNoise(sigma)

X = np.linspace(0, 2, 200)[:, None]
K = cov(X).eval()

plt.plot(MvNormal.dist(mu=np.zeros(200), cov=K).random(size=3).T);
trace 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-829-c9af7d4d5a47> in <module>
      5 K = cov(X).eval()
      6 
----> 7 plt.plot(MvNormal.dist(mu=np.zeros(200), cov=K).random(size=3).T);

<ipython-input-824-191770c5572e> in random(self, point, size)
     99 
    100         mu = broadcast_dist_samples_to(to_shape=output_shape, samples=[mu], size=size)[0]
--> 101         param = np.broadcast_to(param, shape=output_shape + dist_shape[-1:])
    102 
    103         assert mu.shape == output_shape

<__array_function__ internals> in broadcast_to(*args, **kwargs)

~/anaconda3/envs/phd/lib/python3.8/site-packages/numpy/lib/stride_tricks.py in broadcast_to(array, shape, subok)
    178            [1, 2, 3]])
    179     """
--> 180     return _broadcast_to(array, shape, subok=subok, readonly=True)
    181 
    182 

~/anaconda3/envs/phd/lib/python3.8/site-packages/numpy/lib/stride_tricks.py in _broadcast_to(array, shape, subok, readonly)
    121                          'negative')
    122     extras = []
--> 123     it = np.nditer(
    124         (array,), flags=['multi_index', 'refs_ok', 'zerosize_ok'] + extras,
    125         op_flags=['readonly'], itershape=shape, order='C')

ValueError: input operand has more dimensions than allowed by the axis remapping

To make it work:

sigma = 2.0
cov = pm.gp.cov.WhiteNoise(sigma)

X = np.linspace(0, 2, 200)[:, None]
K = cov(X).eval()

plt.plot(MvNormal.dist(mu=np.zeros(200), cov=K, shape=200).random(size=3).T);

I don’t know if this was intentionally defined this way because of some other reason, but it changes the way to use the random method from MvNormal. It is also certainly breaking a lot of examples that use MvNormal. Happy to submit an issue if this is a bug.

Thanks!
Luis

2

Hi @luisroque
When MvNormal is initialized with PyMC3 distributions or symbolic theano tensors as parameters, the shape cannot be inferred. Inferring shape becomes more difficult with operations of symbolic tensors with numpy arrays / theano constants. As a result, providing a shape argument became important to achieve deterministic nature of random method.
The issue has already been reported here asking to raise an error when MvNormal is initialized without shape.

Here, I would certainly ask to open an issue with the list of examples that need to be updated :slight_smile:
Thanks,
Sayam

3

Hi @Sayam753, thank you for the explanation! Just out of curiosity, before this release shape was not required, how was this handled?

I see, thanks for pointing that out! Actually, in my case, the error is different because I’m not supplying observed data. Not sure if you have to grab the error in two places.

Will do! :slight_smile:

1 Like
4

Earlier, we were simply checking if parameters (mu and any one of cov / chol / tau) are drawn with size shape in front. But this is incorrect way of handling shapes because the sample shape (size) can get mixed with distribution shape in case they are equal.

Even if MvNormal is used before 3.10.0 release in model context, there is an error not using shape or observed argument. MvNormal expects at least 1 dimensional tensor for logp, and it receives a scalar if no shape is provided. This leads to an error. So, I think random method has no involvement here.

Thanks

1 Like