Make MvStudentT distribution v4 compatible #4731

Sayam753 · 2021-05-30T12:13:28Z

Thank your for opening a PR!

Depending on what your PR does, here are a few things you might want to address in the description:

Sayam753 · 2021-05-30T13:36:21Z

Looking at the failed pipeline, it seems like it uses cached scipy==1.2.1. But doing a pip install -e . on v4, gives me scipy==1.6.3.

I am using scipy.stats.multivariate_t for drawing random samples which was added in scipy==1.6.0. I think its better to not use a function from a newer version, and revert to using a similar implementation as in master -

https://github.com/pymc-devs/pymc3/blob/d7172c0a1a76301031d1b3b411d00643c416a0c4/pymc3/distributions/multivariate.py#L406-L409

ricardoV94 · 2021-05-30T13:44:29Z

Why not use the newer scipy version?

Sayam753 · 2021-05-30T14:06:23Z

Why not use the newer scipy version?

I think it may lead to changing the requirements file and our pipeline cache. If I already had scipy<1.6.0 cached locally, doing pip install -e . would have installed that earlier version of scipy (which will still be compatible with pymc3 requirements of scipy>=1.2.0)

If using newer scipy versions is feasible with other dependencies' requirements, I am happy to make the required changes :)

pymc3/tests/test_distributions_random.py

pymc3/tests/test_distributions.py

Sayam753 · 2021-06-05T19:56:55Z

scipy.stats.multivariate_t has weird behavior. Its rvs method works only when size parameter is strictly 1 dimensional.

>>> import numpy as np
>>> import scipy
>>> scipy.__version__
'1.6.3'
>>> from scipy.stats import multivariate_normal, multivariate_t
>>> mean, cov = np.zeros(2), np.eye(2)
>>> multivariate_normal.rvs(mean=mean, cov=cov, size=None).shape
(2,)
>>> multivariate_t.rvs(loc=mean, shape=cov, size=None)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 4148, in rvs
    samples = loc + z / np.sqrt(x)[:, None]
IndexError: invalid index to scalar variable.
>>> multivariate_normal.rvs(mean=mean, cov=cov, size=(10, 3)).shape
(10, 3, 2)
>>> multivariate_t.rvs(loc=mean, shape=cov, size=(10, 3))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 4148, in rvs
    samples = loc + z / np.sqrt(x)[:, None]
ValueError: operands could not be broadcast together with shapes (10,3,2) (10,1,3) 
>>> multivariate_t.rvs(loc=mean, shape=cov, size=(3))
array([[ 1.54502738, -0.23866557],
       [-0.79040924, -1.72035808],
       [-0.28060545, -0.69889876]])

This is the reason why many pipelines failed.

Now, I want to try out implementing rng_fn method similar to what we had in v3.

Sayam753 · 2021-06-08T21:11:20Z

I am on latest main branch@ 4de95c3 and trying to understand how new API for shapes work.

A small snippet -

>>> import pymc3 as pm
You are running the v4 development version of PyMC3 which currently still lacks key features. You probably want to use the stable v3 instead which you can either install via conda or find on the v3 GitHub branch: https://github.com/pymc-devs/pymc3/tree/v3
>>> import numpy as np
>>> dist = pm.MvNormal.dist(mu=np.zeros(1), cov=np.eye(3), shape=3)
>>> dist.shape.eval()
array([1])

I am finding it hard to understand why the shape is [1]
Shouldn't it be [3]?

And some fancy errors during dist.eval()

Traceback

>>> dist.eval()
Traceback (most recent call last):
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/compile/function/types.py", line 977, in __call__
    if output_subset is None
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/graph/op.py", line 495, in rval
    r = p(n, [x[0] for x in i], o)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/random/op.py", line 413, in perform
    smpl_val = self.rng_fn(rng, *(args + [size]))
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/random/basic.py", line 288, in rng_fn
    return safe_multivariate_normal(mean, cov, size=size, rng=rng)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/random/basic.py", line 239, in safe_multivariate_normal
    stats.multivariate_normal(mean=mean, cov=cov, allow_singular=True).rvs(
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 369, in __call__
    seed=seed)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 744, in __init__
    None, mean, cov)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 409, in _process_parameters
    cov.shape = (1, 1)
ValueError: cannot reshape array of size 9 into shape (1,1)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/graph/basic.py", line 555, in eval
    rval = self._fn_cache[inputs](*args)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/compile/function/types.py", line 993, in __call__
    storage_map=getattr(self.fn, "storage_map", None),
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/link/utils.py", line 521, in raise_with_op
    raise exc_value.with_traceback(exc_trace)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/compile/function/types.py", line 977, in __call__
    if output_subset is None
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/graph/op.py", line 495, in rval
    r = p(n, [x[0] for x in i], o)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/random/op.py", line 413, in perform
    smpl_val = self.rng_fn(rng, *(args + [size]))
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/random/basic.py", line 288, in rng_fn
    return safe_multivariate_normal(mean, cov, size=size, rng=rng)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/random/basic.py", line 239, in safe_multivariate_normal
    stats.multivariate_normal(mean=mean, cov=cov, allow_singular=True).rvs(
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 369, in __call__
    seed=seed)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 744, in __init__
    None, mean, cov)
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/scipy/stats/_multivariate.py", line 409, in _process_parameters
    cov.shape = (1, 1)
ValueError: cannot reshape array of size 9 into shape (1,1)
Apply node that caused the error: multivariate_normal_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FB57E736380>), TensorConstant{[]}, TensorConstant{11}, TensorConstant{(1,) of 0.0}, TensorConstant{[[1. 0. 0...0. 0. 1.]]})
Toposort index: 0
Inputs types: [RandomStateType, TensorType(int64, vector), TensorType(int64, scalar), TensorType(float64, (True,)), TensorType(float64, matrix)]
Inputs shapes: ['No shapes', (0,), (), (1,), (3, 3)]
Inputs strides: ['No strides', (8,), (), (8,), (24, 8)]
Inputs values: [RandomState(MT19937) at 0x7FB57E736380, array([], dtype=int64), array(11), array([0.]), 'not shown']
Outputs clients: [[], ['output']]

Backtrace when the node is created (use Aesara flag traceback__limit=N to make it longer):
  File "<stdin>", line 1, in <module>
  File "/Users/sayam/Desktop/pymc/pymc3/pymc3/distributions/multivariate.py", line 226, in dist
    return super().dist([mu, cov], **kwargs)
  File "/Users/sayam/Desktop/pymc/pymc3/pymc3/distributions/distribution.py", line 285, in dist
    rv_out = cls.rv_op(*dist_params, size=create_size, **kwargs)

HINT: Use the Aesara flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

cc @michaelosthege

michaelosthege · 2021-06-08T21:25:25Z

>>> dist = pm.MvNormal.dist(mu=np.zeros(1), cov=np.eye(3), shape=3)
>>> dist.shape.eval()
array([1])
I am finding it hard to understand why the shape is [1]
Shouldn't it be [3]?

This parametrization is incoherent. However, this is a scenario that motivates automatically adding a SpecifyShape Op which I had implemented but then took out, because I couldn't come up with a test case to trigger it.

Your shape kwarg has no effect, because the last aesara.tensor.random.basic.multivariate_normal.ndim_supp dimensions are shaved off.
Internally it becomes:

>>> import aesara
>>> aesara.__version__
'2.0.8'
>>> import aesara.tensor as at
>>> import aesara.tensor.random.basic as ar
>>> import numpy as np
>>> ar.multivariate_normal(mean=np.zeros(1), cov=np.eye(3)).shape.eval()
array([1], dtype=int64)

The mean/mu parameter to the MvNormal should include the support dimensions. So for example:

>>> ar.multivariate_normal(mean=[0,0,0], cov=np.eye(3)).shape.eval()
array([3], dtype=int64)

One could argue that the Aesara implementation shouldn't have downcasted cov just because mean had shape (1,). At least not without printing a warning.
cc @brandonwillard

Sayam753 · 2021-06-08T22:06:02Z

The mean/mu parameter to the MvNormal should include the support dimensions

I guess this is due to rep_param_idx parameter used to infer distribution shape in aesara.

One could argue that the Aesara implementation shouldn't have downcasted cov just because mean had shape (1,). At least not without printing a warning.

How hard would it be to loop over all distribution parameters to infer distribution shape, instead of relying on one parameter?

brandonwillard · 2021-06-09T02:11:20Z

However, this is a scenario that motivates automatically adding a SpecifyShape Op which I had implemented but then took out, because I couldn't come up with a test case to trigger it.

SpecifyShape wouldn't help at all.

The mean/mu parameter to the MvNormal should include the support dimensions. So for example:
>>> ar.multivariate_normal(mean=[0,0,0], cov=np.eye(3)).shape.eval()
array([3], dtype=int64)

This has nothing to do with "including the support dimensions"; the actual issue is that the parameters specify conflicting lengths/sizes across the relevant dimensions. Your follow-up example isn't including anything extra, it's only fixing the conflicting dimensions by making the length of the mean match the lengths of the covariance dimensions.

Consider what [0, 0, 0] actually is: a one-dimensional vector of length three. It's not a three-dimensional array.

I think a lot of confusion regarding shape, size, etc., might be stemming from this fundamental misconception regarding dimensions and their sizes/lengths.

One could argue that the Aesara implementation shouldn't have downcasted cov just because mean had shape (1,). At least not without printing a warning.

What "downcasting"? The dtype looks correct.

The result itself appears to be due to a missing check in the shape inference. That's all.

If you want to understand why, start by computing the result itself, and not just its shape. You'll see the underlying error that the shape inference logic missed:

>>> ar.multivariate_normal(mean=np.zeros(1), cov=np.eye(3)).eval()
...
ValueError: cannot reshape array of size 9 into shape (1,1)

This is the same error that appears in the original problem. In other words, the shape parameter isn't even relevant, nor is the MvNormal wrapper.

As usual, it's always good to see what the underlying NumPy function would do, since these Ops are modeling/wrapping those functions:

>>> np.random.multivariate_normal(mean=np.zeros(1), cov=np.eye(3))
ValueError: mean and cov must have same length

As expected, NumPy fails as well; however, it gives a much better error message.

@Sayam753, feel free to open an issue in Aesara to add that missing check to the shape inference steps; otherwise, your example specifies a mathematically invalid distribution, because the mean and covariance are incompatible. Use mu=np.zeros(3) if you want a multivariate normal that produces vectors of length three instead.

michaelosthege · 2021-06-10T09:43:21Z

SpecifyShape wouldn't help at all.

IIRC I once had a version where SpecifyShape worked and made sense (when working on #4552). I don't remeber the details or what changed though.
Brandon is right: The ValueError: operands could not be broadcast together is raised before a SpecifyShape can get involved.

What "downcasting"? The dtype looks correct.

Sorry, I mean down_sizing_. (The bug in the shape inference.)

This has nothing to do with "including the support dimensions"; the actual issue is that the parameters specify conflicting lengths/sizes across the relevant dimensions.

Please excuse that my language wasn't exactly on point. The phrase "including the support dimensions" is more from the perspective of wanting to imply more dimensions from the RV parameters, which wasn't the intention in the first place.

Anyways. @Sayam753 let us know if you need assistance with the remaining ToDo for this PR.

Sayam753 · 2021-06-11T19:33:08Z

Thanks @michaelosthege and @brandonwillard for clearing the doubts above.

A little interesting behaviour of MvNormal.logp as it works even when mean and cov are not having the same length -

>>> import numpy as np
>>> import pymc3 as pm
You are running the v4 development version of PyMC3 which currently still lacks key features. You probably want to use the stable v3 instead which you can either install via conda or find on the v3 GitHub branch: https://github.com/pymc-devs/pymc3/tree/v3
>>> mu, cov = np.ones(1), np.eye(3)
>>> pm.MvNormal.logp(value=np.random.randn(20, 3), mu=mu, cov=cov).eval()
array([ -3.67340545,  -7.19130309,  -3.23923742,  -5.15713797,
       -17.29246223, -10.32441825,  -4.2939353 ,  -5.66928884,
        -4.05121295,  -4.58167586,  -5.50445437,  -4.20362596,
        -6.12579132,  -8.63937152,  -4.656641  ,  -9.63474273,
        -4.35076017,  -3.33788195,  -3.50893476,  -4.77016247])

ricardoV94 · 2021-06-11T22:49:36Z

Calling the logp method directly does not go through the random variable op, so no shape or size checks would be done (it's just a normal aesara function call). As long as the function broadcasts well you wouldn't see any error.

Sayam753 · 2021-06-12T18:26:25Z

I am trying to understand how logpt works under the hood. This will help in writing a better solution for the test case pymc3/tests/test_distributions.py::TestBugfixes::test_issue_3051

My question is - why cannot I use filter_variable method to broadcast a vector of shape (3) and a matrix of shape (20, 3) in aesara?

>>> import numpy as np
>>> import aesara.tensor as at
>>> import pymc3 as pm
You are running the v4 development version of PyMC3 which currently still lacks key features. You probably want to use the stable v3 instead which you can either install via conda or find on the v3 GitHub branch: https://github.com/pymc-devs/pymc3/tree/v3
>>> mu, cov = np.ones(3), np.eye(3)
>>> dist = pm.MvNormal.dist(mu=mu, cov=cov)
>>> X = np.random.normal(size=(20, 3))

>>> X = at.as_tensor(X)
>>> dist.type, dist.dtype, dist.ndim, dist.shape.eval()
(TensorType(float64, vector), 'float64', 1, array([3]))
>>> X.type, X.dtype, X.ndim, X.shape.eval()
(TensorType(float64, matrix), 'float64', 2, array([20,  3]))

>>> dist.type.filter_variable(X)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/sayam/miniconda3/envs/pymc_dev/lib/python3.7/site-packages/aesara/tensor/type.py", line 261, in filter_variable
    f"Cannot convert Type {other.type} "
TypeError: Cannot convert Type TensorType(float64, matrix) (of Variable TensorConstant{[[ 0.30673..15832267]]}) into Type TensorType(float64, vector). You can try to manually convert TensorConstant{[[ 0.30673..15832267]]} into a TensorType(float64, vector).

cc @brandonwillard

Sayam753 · 2021-06-13T03:20:02Z

The pipelines are green. So, the PR is ready for another round of review.

pymc3/distributions/multivariate.py

ricardoV94 · 2021-06-13T07:42:42Z

In the meantime, were there no other tests marked as xfail because they relied on this distribution?

Sayam753 · 2021-06-13T13:19:35Z

I think all the MvStudentT related tests in test_distributions.py and test_distributions_random.py have been refactored/added in this PR.

Thinking more about the test cases, MvStudentTRandomWalk's rng method must depend on MvStudentT rng method. AFAIK the RandomWalk distributions have not been refactored yet and including them here will be an overkill for this PR.

pymc3/distributions/multivariate.py

ricardoV94

I left some minor suggestions. Otherwise it looks good to me. Great work :)

Co-authored-by: Ricardo <[email protected]>

ricardoV94

LGTM!

ricardoV94 · 2021-06-14T15:44:02Z

Thanks @Sayam753, great work.

Sayam753 · 2021-06-14T15:54:52Z

Thanks @ricardoV94 for a thorough review. I do want to follow up for a better fix for pymc3/tests/test_distributions.py::TestBugfixes::test_issue_3051 and will open a related issue :)

ricardoV94 reviewed May 30, 2021

View reviewed changes

pymc3/tests/test_distributions_random.py Show resolved Hide resolved

ricardoV94 reviewed May 30, 2021

View reviewed changes

pymc3/tests/test_distributions.py Outdated Show resolved Hide resolved

ricardoV94 mentioned this pull request May 30, 2021

Port remaining distributions to v4 #4686

Closed

26 tasks

Sayam753 changed the base branch from v4 to main June 5, 2021 15:12

Sayam753 force-pushed the add_mvstudentt branch from 6239763 to 61a5b8f Compare June 5, 2021 15:13

Sayam753 mentioned this pull request Jun 6, 2021

Incorrect shape handling in scipy.stat.multivariate_t.rvs method scipy/scipy#14189

Closed

ricardoV94 requested changes Jun 13, 2021

View reviewed changes

pymc3/distributions/multivariate.py Outdated Show resolved Hide resolved

ricardoV94 reviewed Jun 13, 2021

View reviewed changes

pymc3/distributions/multivariate.py Outdated Show resolved Hide resolved

ricardoV94 reviewed Jun 13, 2021

View reviewed changes

pymc3/distributions/multivariate.py Outdated Show resolved Hide resolved

ricardoV94 reviewed Jun 13, 2021

View reviewed changes

Sayam753 added 5 commits June 13, 2021 20:53

Make MvStudentT distribution v4 compatible

e1d9c69

Enable pymc3/tests/test_distributions.py::TestBugfixes::test_issue_3051

1425508

refactor tests

8eb5623

simplify rng_fn method and refactor tests

a553c55

Refactor rng method to compute sqrt over chi2_samples, one step earlier

91e2dd1

Sayam753 and others added 3 commits June 13, 2021 20:53

Modify rng method to use aesara multivariate_normal.rng_fn

eb3901d

Add suggestions

c98b214

Co-authored-by: Ricardo <[email protected]>

Add a check against negative support of degrees of freedom

de30b28

Co-authored-by: Ricardo <[email protected]>

Sayam753 force-pushed the add_mvstudentt branch from 5ab17d0 to de30b28 Compare June 13, 2021 15:25

ricardoV94 approved these changes Jun 13, 2021

View reviewed changes

ricardoV94 merged commit f8a5b91 into pymc-devs:main Jun 14, 2021

ricardoV94 mentioned this pull request Jun 18, 2021

[WIP] Porting kroneckernormal distribution to v4 #4774

Merged

Sayam753 deleted the add_mvstudentt branch August 16, 2021 13:25

kc611 mentioned this pull request Sep 22, 2021

Refactored LKJCorr and _LKJCholeskyCov #4784

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Make MvStudentT distribution v4 compatible #4731

Make MvStudentT distribution v4 compatible #4731

Sayam753 commented May 30, 2021 •

edited

Loading

Sayam753 commented May 30, 2021

ricardoV94 commented May 30, 2021

Sayam753 commented May 30, 2021

Sayam753 commented Jun 5, 2021

Sayam753 commented Jun 8, 2021 •

edited

Loading

michaelosthege commented Jun 8, 2021

Sayam753 commented Jun 8, 2021 •

edited

Loading

brandonwillard commented Jun 9, 2021 •

edited

Loading

michaelosthege commented Jun 10, 2021

Sayam753 commented Jun 11, 2021

ricardoV94 commented Jun 11, 2021 •

edited

Loading

Sayam753 commented Jun 12, 2021

Sayam753 commented Jun 13, 2021

ricardoV94 commented Jun 13, 2021 •

edited

Loading

Sayam753 commented Jun 13, 2021 •

edited

Loading

ricardoV94 left a comment •

edited

Loading

ricardoV94 left a comment

ricardoV94 commented Jun 14, 2021

Sayam753 commented Jun 14, 2021

Make MvStudentT distribution v4 compatible #4731

Make MvStudentT distribution v4 compatible #4731

Conversation

Sayam753 commented May 30, 2021 • edited Loading

Sayam753 commented May 30, 2021

ricardoV94 commented May 30, 2021

Sayam753 commented May 30, 2021

Sayam753 commented Jun 5, 2021

Sayam753 commented Jun 8, 2021 • edited Loading

michaelosthege commented Jun 8, 2021

Sayam753 commented Jun 8, 2021 • edited Loading

brandonwillard commented Jun 9, 2021 • edited Loading

michaelosthege commented Jun 10, 2021

Sayam753 commented Jun 11, 2021

ricardoV94 commented Jun 11, 2021 • edited Loading

Sayam753 commented Jun 12, 2021

Sayam753 commented Jun 13, 2021

ricardoV94 commented Jun 13, 2021 • edited Loading

Sayam753 commented Jun 13, 2021 • edited Loading

ricardoV94 left a comment • edited Loading

Choose a reason for hiding this comment

ricardoV94 left a comment

Choose a reason for hiding this comment

ricardoV94 commented Jun 14, 2021

Sayam753 commented Jun 14, 2021

Sayam753 commented May 30, 2021 •

edited

Loading

Sayam753 commented Jun 8, 2021 •

edited

Loading

Sayam753 commented Jun 8, 2021 •

edited

Loading

brandonwillard commented Jun 9, 2021 •

edited

Loading

ricardoV94 commented Jun 11, 2021 •

edited

Loading

ricardoV94 commented Jun 13, 2021 •

edited

Loading

Sayam753 commented Jun 13, 2021 •

edited

Loading

ricardoV94 left a comment •

edited

Loading