From fe1128b9ffe31355eaf13ccd5fa18939af3463dd Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Thu, 11 Nov 2021 22:00:31 +0530 Subject: [PATCH 01/11] Add MvStudentT moment --- pymc/distributions/multivariate.py | 10 +++++++++- pymc/tests/test_distributions_moments.py | 16 ++++++++++++++++ 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index cfa372e037..7b21a8dda3 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -46,7 +46,7 @@ from pymc.distributions.continuous import ChiSquared, Normal, assert_negative_support from pymc.distributions.dist_math import bound, factln, logpow, multigammaln from pymc.distributions.distribution import Continuous, Discrete -from pymc.distributions.shape_utils import broadcast_dist_samples_to, to_tuple +from pymc.distributions.shape_utils import broadcast_dist_samples_to, rv_size_is_none, to_tuple from pymc.math import kron_diag, kron_dot __all__ = [ @@ -343,6 +343,14 @@ def dist(cls, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None, lower=True assert_negative_support(nu, "nu", "MvStudentT") return super().dist([nu, mu, cov], **kwargs) + def get_moment(rv, size, nu, mu, cov): + moment = mu + if not rv_size_is_none(size): + if isinstance(size, int): + size = (size,) + moment = at.full((*size, mu.size), moment) + return moment + def logp(value, nu, mu, cov): """ Calculate log-probability of Multivariate Student's T distribution diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index f52ca8baa0..4b5a6ca92d 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -25,6 +25,7 @@ Laplace, Logistic, LogNormal, + MvStudentT, NegativeBinomial, Normal, Pareto, @@ -612,3 +613,18 @@ def test_discrete_uniform_moment(lower, upper, size, expected): with Model() as model: DiscreteUniform("x", lower=lower, upper=upper, size=size) assert_moment_is_expected(model, expected) + + +@pytest.mark.parametrize( + "nu, mu, cov, size, expected", + [ + (2, np.ones(1), np.eye(1), None, np.ones(1)), + (2, np.ones(2), np.eye(2), None, np.ones(2)), + (2, np.ones(2), np.eye(2), 2, np.full((2, 2), np.ones(2))), + (1, np.ones(2), np.eye(2), (2, 5), np.full((2, 5, 2), np.ones(2))), + ], +) +def test_mvstudentt_moment(nu, mu, cov, size, expected): + with Model() as model: + MvStudentT("x", nu=nu, mu=mu, cov=cov, size=size) + assert_moment_is_expected(model, expected) From 6ea3b41c2a405745bad8853f46de36449a1c4097 Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Thu, 11 Nov 2021 22:19:15 +0530 Subject: [PATCH 02/11] fix pre-commit --- pymc/distributions/multivariate.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index 7b21a8dda3..c8599b814e 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -46,7 +46,11 @@ from pymc.distributions.continuous import ChiSquared, Normal, assert_negative_support from pymc.distributions.dist_math import bound, factln, logpow, multigammaln from pymc.distributions.distribution import Continuous, Discrete -from pymc.distributions.shape_utils import broadcast_dist_samples_to, rv_size_is_none, to_tuple +from pymc.distributions.shape_utils import ( + broadcast_dist_samples_to, + rv_size_is_none, + to_tuple, +) from pymc.math import kron_diag, kron_dot __all__ = [ From 7b3e63328fb9361b85f8e8c8f6d74a96efa2e4d0 Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Fri, 12 Nov 2021 08:47:14 +0530 Subject: [PATCH 03/11] Add MatrixNormal moment --- pymc/distributions/multivariate.py | 8 ++++++++ pymc/tests/test_distributions_moments.py | 26 +++++++++++++++++++++--- 2 files changed, 31 insertions(+), 3 deletions(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index c8599b814e..40c8e20058 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -1719,6 +1719,14 @@ def dist( return super().dist([mu, rowchol_cov, colchol_cov], **kwargs) + def get_moment(rv, size, mu, rowchol, colchol): + moment = mu + if not rv_size_is_none(size): + if isinstance(size, int): + size = (size,) + moment = at.full((*size, mu.size), moment) + return moment + def logp(value, mu, rowchol, colchol): """ Calculate log-probability of Matrix-valued Normal distribution diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index 4b5a6ca92d..351e6863ce 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -25,6 +25,7 @@ Laplace, Logistic, LogNormal, + MatrixNormal, MvStudentT, NegativeBinomial, Normal, @@ -615,16 +616,35 @@ def test_discrete_uniform_moment(lower, upper, size, expected): assert_moment_is_expected(model, expected) +rand1d = np.random.rand(2) + + @pytest.mark.parametrize( "nu, mu, cov, size, expected", [ (2, np.ones(1), np.eye(1), None, np.ones(1)), - (2, np.ones(2), np.eye(2), None, np.ones(2)), - (2, np.ones(2), np.eye(2), 2, np.full((2, 2), np.ones(2))), - (1, np.ones(2), np.eye(2), (2, 5), np.full((2, 5, 2), np.ones(2))), + (2, rand1d, np.eye(2), None, rand1d), + (2, rand1d, np.eye(2), 2, np.full((2, 2), rand1d)), + (1, rand1d, np.eye(2), (2, 5), np.full((2, 5, 2), rand1d)), ], ) def test_mvstudentt_moment(nu, mu, cov, size, expected): with Model() as model: MvStudentT("x", nu=nu, mu=mu, cov=cov, size=size) assert_moment_is_expected(model, expected) + + +rand2d = np.random.rand(2, 3) + + +@pytest.mark.parametrize( + "mu, rowchol, colchol, expected", + [ + (np.ones((1, 1)), np.eye(1), np.eye(1), np.ones((1, 1))), + (rand2d, np.eye(2), np.eye(3), rand2d), + ], +) +def test_matrixnormal_moment(mu, rowchol, colchol, expected): + with Model() as model: + MatrixNormal("x", mu=mu, rowchol=rowchol, colchol=colchol) + assert_moment_is_expected(model, expected) From 97d0ce542ec27e91e686f31302196d073ae86563 Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Sun, 14 Nov 2021 14:55:12 +0530 Subject: [PATCH 04/11] fix moments and cover with tests --- pymc/distributions/multivariate.py | 16 +++++++--------- pymc/tests/test_distributions_moments.py | 21 ++++++++++++--------- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index 40c8e20058..e66f64ba41 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -350,9 +350,8 @@ def dist(cls, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None, lower=True def get_moment(rv, size, nu, mu, cov): moment = mu if not rv_size_is_none(size): - if isinstance(size, int): - size = (size,) - moment = at.full((*size, mu.size), moment) + moment_size = at.concatenate([size, mu.shape]) + moment = at.full(moment_size, moment) return moment def logp(value, nu, mu, cov): @@ -1668,8 +1667,8 @@ def dist( cholesky = Cholesky(lower=True, on_error="raise") - if kwargs.get("size", None) is not None: - raise NotImplementedError("MatrixNormal doesn't support size argument") + # if kwargs.get("size", None) is not None: + # raise NotImplementedError("MatrixNormal doesn't support size argument") if "shape" in kwargs: kwargs.pop("shape") @@ -1720,11 +1719,10 @@ def dist( return super().dist([mu, rowchol_cov, colchol_cov], **kwargs) def get_moment(rv, size, mu, rowchol, colchol): - moment = mu + output_shape = (rowchol.shape[0], colchol.shape[0]) if not rv_size_is_none(size): - if isinstance(size, int): - size = (size,) - moment = at.full((*size, mu.size), moment) + output_shape = at.concatenate([size, output_shape]) + moment = at.full(output_shape, mu) return moment def logp(value, mu, rowchol, colchol): diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index 351e6863ce..d4284c1bfb 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -617,6 +617,7 @@ def test_discrete_uniform_moment(lower, upper, size, expected): rand1d = np.random.rand(2) +rand2d = np.random.rand(2, 3) @pytest.mark.parametrize( @@ -625,7 +626,9 @@ def test_discrete_uniform_moment(lower, upper, size, expected): (2, np.ones(1), np.eye(1), None, np.ones(1)), (2, rand1d, np.eye(2), None, rand1d), (2, rand1d, np.eye(2), 2, np.full((2, 2), rand1d)), - (1, rand1d, np.eye(2), (2, 5), np.full((2, 5, 2), rand1d)), + (2, rand2d, np.eye(3), 2, np.full((2, 2, 3), rand2d)), + (2, rand1d, np.eye(2), (2, 5), np.full((2, 5, 2), rand1d)), + (2, rand2d, np.eye(3), (2, 5), np.full((2, 5, 2, 3), rand2d)), ], ) def test_mvstudentt_moment(nu, mu, cov, size, expected): @@ -634,17 +637,17 @@ def test_mvstudentt_moment(nu, mu, cov, size, expected): assert_moment_is_expected(model, expected) -rand2d = np.random.rand(2, 3) - - @pytest.mark.parametrize( - "mu, rowchol, colchol, expected", + "mu, rowchol, colchol, size, expected", [ - (np.ones((1, 1)), np.eye(1), np.eye(1), np.ones((1, 1))), - (rand2d, np.eye(2), np.eye(3), rand2d), + (np.ones((1, 1)), np.eye(1), np.eye(1), None, np.ones((1, 1))), + (np.ones((1, 1)), np.eye(2), np.eye(3), None, np.ones((2, 3))), + (rand2d, np.eye(2), np.eye(3), None, rand2d), + (rand2d, np.eye(2), np.eye(3), 2, np.full((2, 2, 3), rand2d)), + (rand2d, np.eye(2), np.eye(3), (2, 5), np.full((2, 5, 2, 3), rand2d)), ], ) -def test_matrixnormal_moment(mu, rowchol, colchol, expected): +def test_matrixnormal_moment(mu, rowchol, colchol, size, expected): with Model() as model: - MatrixNormal("x", mu=mu, rowchol=rowchol, colchol=colchol) + MatrixNormal("x", mu=mu, rowchol=rowchol, colchol=colchol, size=size) assert_moment_is_expected(model, expected) From 964c6b493df8af2ede9510b181fd8f14fc0e7cf3 Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Sun, 14 Nov 2021 19:46:23 +0530 Subject: [PATCH 05/11] disallow size in MatrixNormal --- pymc/distributions/multivariate.py | 4 ++-- pymc/tests/test_distributions_moments.py | 14 +++++++++++--- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index e66f64ba41..e38f19cb6f 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -1667,8 +1667,8 @@ def dist( cholesky = Cholesky(lower=True, on_error="raise") - # if kwargs.get("size", None) is not None: - # raise NotImplementedError("MatrixNormal doesn't support size argument") + if kwargs.get("size", None) is not None: + raise NotImplementedError("MatrixNormal doesn't support size argument") if "shape" in kwargs: kwargs.pop("shape") diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index d4284c1bfb..6296ce5e46 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -637,6 +637,12 @@ def test_mvstudentt_moment(nu, mu, cov, size, expected): assert_moment_is_expected(model, expected) +def check_matrixnormal_moment(mu, rowchol, colchol, size, expected): + with Model() as model: + MatrixNormal("x", mu=mu, rowchol=rowchol, colchol=colchol, size=size) + assert_moment_is_expected(model, expected) + + @pytest.mark.parametrize( "mu, rowchol, colchol, size, expected", [ @@ -648,6 +654,8 @@ def test_mvstudentt_moment(nu, mu, cov, size, expected): ], ) def test_matrixnormal_moment(mu, rowchol, colchol, size, expected): - with Model() as model: - MatrixNormal("x", mu=mu, rowchol=rowchol, colchol=colchol, size=size) - assert_moment_is_expected(model, expected) + if size is None: + check_matrixnormal_moment(mu, rowchol, colchol, size, expected) + else: + with pytest.raises(NotImplementedError): + check_matrixnormal_moment(mu, rowchol, colchol, size, expected) From 59f44d1f20af35f44100ce7632c9b8b2b0856bad Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Sun, 14 Nov 2021 21:21:03 +0530 Subject: [PATCH 06/11] fix pre-commit --- pymc/tests/test_distributions_moments.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index 91fb1f4756..e4df6b17db 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -31,8 +31,8 @@ LogitNormal, LogNormal, MatrixNormal, - MvStudentT, Moyal, + MvStudentT, NegativeBinomial, Normal, Pareto, From a26560d01edddf668232c1862fe2aed26c7b3316 Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Sun, 14 Nov 2021 22:14:55 +0530 Subject: [PATCH 07/11] empty commit to rerun the tests From d2537ed514a2afd1ce0bc4e8ee713f7bec435ced Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Sun, 14 Nov 2021 22:21:27 +0530 Subject: [PATCH 08/11] empty commit to rerun the tests From 3e0c66e5ead9d600f23e690655f13a6883a7645c Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Wed, 17 Nov 2021 12:36:54 +0530 Subject: [PATCH 09/11] Fix `nu` broadcasting --- pymc/distributions/multivariate.py | 9 +++++---- pymc/tests/test_distributions_moments.py | 1 + 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index b3cc7f39f8..06ecaeb05e 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -348,10 +348,11 @@ def dist(cls, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None, lower=True return super().dist([nu, mu, cov], **kwargs) def get_moment(rv, size, nu, mu, cov): - moment = mu - if not rv_size_is_none(size): - moment_size = at.concatenate([size, mu.shape]) - moment = at.full(moment_size, moment) + if rv_size_is_none(size): + moment_size = at.concatenate([nu.shape, mu.shape]) + else: + moment_size = at.concatenate([size, nu.shape, mu.shape]) + moment = at.full(moment_size, mu) return moment def logp(value, nu, mu, cov): diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index e4df6b17db..2e4c1c1126 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -781,6 +781,7 @@ def test_moyal_moment(mu, sigma, size, expected): (2, rand2d, np.eye(3), 2, np.full((2, 2, 3), rand2d)), (2, rand1d, np.eye(2), (2, 5), np.full((2, 5, 2), rand1d)), (2, rand2d, np.eye(3), (2, 5), np.full((2, 5, 2, 3), rand2d)), + (np.array([3, 4]), np.ones(2), np.eye(2), None, np.ones((2, 2))), ], ) def test_mvstudentt_moment(nu, mu, cov, size, expected): From 766a25d047cf4908bcfb5517e7219bae71828aaf Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Wed, 17 Nov 2021 12:53:52 +0530 Subject: [PATCH 10/11] Fix pre-commit --- pymc/tests/test_distributions_moments.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index 32cfafeac5..e1687b1501 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -806,7 +806,6 @@ def test_moyal_moment(mu, sigma, size, expected): assert_moment_is_expected(model, expected) - rand1d = np.random.rand(2) rand2d = np.random.rand(2, 3) From e23a67281feeedbf559bf9f27e5dc0821be0051d Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Sat, 20 Nov 2021 22:16:32 +0530 Subject: [PATCH 11/11] disable broadcasting --- pymc/distributions/multivariate.py | 11 +++++------ pymc/tests/test_distributions_moments.py | 4 ++-- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/pymc/distributions/multivariate.py b/pymc/distributions/multivariate.py index cd1c114e89..4dbc1ea9b8 100644 --- a/pymc/distributions/multivariate.py +++ b/pymc/distributions/multivariate.py @@ -355,11 +355,10 @@ def dist(cls, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None, lower=True return super().dist([nu, mu, cov], **kwargs) def get_moment(rv, size, nu, mu, cov): - if rv_size_is_none(size): - moment_size = at.concatenate([nu.shape, mu.shape]) - else: - moment_size = at.concatenate([size, nu.shape, mu.shape]) - moment = at.full(moment_size, mu) + moment = mu + if not rv_size_is_none(size): + moment_size = at.concatenate([size, moment.shape]) + moment = at.full(moment_size, moment) return moment def logp(value, nu, mu, cov): @@ -700,7 +699,7 @@ def dist(cls, eta, cutpoints, n, *args, **kwargs): class OrderedMultinomial: - R""" + r""" Wrapper class for Ordered Multinomial distributions. Useful for regression on ordinal data whose values range diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py index 5325560824..8b851c0044 100644 --- a/pymc/tests/test_distributions_moments.py +++ b/pymc/tests/test_distributions_moments.py @@ -842,10 +842,10 @@ def test_moyal_moment(mu, sigma, size, expected): (2, np.ones(1), np.eye(1), None, np.ones(1)), (2, rand1d, np.eye(2), None, rand1d), (2, rand1d, np.eye(2), 2, np.full((2, 2), rand1d)), - (2, rand2d, np.eye(3), 2, np.full((2, 2, 3), rand2d)), (2, rand1d, np.eye(2), (2, 5), np.full((2, 5, 2), rand1d)), + (2, rand2d, np.eye(3), None, rand2d), + (2, rand2d, np.eye(3), 2, np.full((2, 2, 3), rand2d)), (2, rand2d, np.eye(3), (2, 5), np.full((2, 5, 2, 3), rand2d)), - (np.array([3, 4]), np.ones(2), np.eye(2), None, np.ones((2, 2))), ], ) def test_mvstudentt_moment(nu, mu, cov, size, expected):