From 8e6f1316020c8e99a2613abc96aac8392cdd1196 Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Thu, 6 Jun 2019 10:40:18 +0200 Subject: [PATCH 1/3] Changed Triangular, TruncatedNormal, Rice and ZeroInflatedNegativeBinomial random method. The math operations between the distribution parameters was moved to _random, and all parameter broadcasting is handled in generate_samples. Added tests. --- pymc3/distributions/continuous.py | 70 ++++++++-- pymc3/distributions/discrete.py | 22 +++- pymc3/tests/test_distributions_random.py | 160 +++++++++++++++++++++++ 3 files changed, 235 insertions(+), 17 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 96a5b554905..343598be846 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -664,16 +664,34 @@ def random(self, point=None, size=None): ------- array """ - mu_v, std_v, a_v, b_v = draw_values( - [self.mu, self.sigma, self.lower, self.upper], point=point, size=size) - return generate_samples(stats.truncnorm.rvs, - a=(a_v - mu_v)/std_v, - b=(b_v - mu_v) / std_v, - loc=mu_v, - scale=std_v, - dist_shape=self.shape, - size=size, - ) + mu, sigma, lower, upper = draw_values( + [self.mu, self.sigma, self.lower, self.upper], + point=point, + size=size + ) + return generate_samples( + self._random, + mu=mu, + sigma=sigma, + lower=lower, + upper=upper, + dist_shape=self.shape, + size=size, + ) + + def _random(self, mu, sigma, lower, upper, size): + """ Wrapper around stats.truncnorm.rvs that converts TruncatedNormal's + parametrization to scipy.truncnorm. All parameter arrays should have + been broadcasted properly by generate_samples at this point and size is + the scipy.rvs representation. + """ + return stats.truncnorm.rvs( + a=(lower - mu) / sigma, + b=(upper - mu) / sigma, + loc=mu, + scale=sigma, + size=size, + ) def logp(self, value): """ @@ -3582,11 +3600,23 @@ def random(self, point=None, size=None): """ c, lower, upper = draw_values([self.c, self.lower, self.upper], point=point, size=size) - scale = upper - lower - c_ = (c - lower) / scale - return generate_samples(stats.triang.rvs, c=c_, loc=lower, scale=scale, + return generate_samples(self._random, c=c, lower=lower, upper=upper, size=size, dist_shape=self.shape) + def _random(self, c, lower, upper, size): + """ Wrapper around stats.triang.rvs that converts Triangular's + parametrization to scipy.triang. All parameter arrays should have + been broadcasted properly by generate_samples at this point and size is + the scipy.rvs representation. + """ + scale = upper - lower + return stats.triang.rvs( + c=(c - lower) / scale, + loc=lower, + scale=scale, + size=size, + ) + def logp(self, value): """ Calculate log-probability of Triangular distribution at specified value. @@ -3875,9 +3905,21 @@ def random(self, point=None, size=None): """ nu, sigma = draw_values([self.nu, self.sigma], point=point, size=size) - return generate_samples(stats.rice.rvs, b=nu / sigma, scale=sigma, loc=0, + return generate_samples(self._random, nu=nu, sigma=sigma, dist_shape=self.shape, size=size) + def _random(self, nu, sigma, size): + """ Wrapper around stats.rice.rvs that converts Rice's + parametrization to scipy.rice. All parameter arrays should have + been broadcasted properly by generate_samples at this point and size is + the scipy.rvs representation. + """ + return stats.rice.rvs( + b=nu / sigma, + scale=sigma, + size=size, + ) + def logp(self, value): """ Calculate log-probability of Rice distribution at specified value. diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index dc40745fab6..8fae8b8cceb 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -1427,13 +1427,29 @@ def random(self, point=None, size=None): """ mu, alpha, psi = draw_values( [self.mu, self.alpha, self.psi], point=point, size=size) - g = generate_samples(stats.gamma.rvs, alpha, scale=mu / alpha, - dist_shape=self.shape, - size=size) + g = generate_samples( + self._random, + mu=mu, + alpha=alpha, + dist_shape=self.shape, + size=size + ) g[g == 0] = np.finfo(float).eps # Just in case g, psi = broadcast_distribution_samples([g, psi], size=size) return stats.poisson.rvs(g) * (np.random.random(g.shape) < psi) + def _random(self, mu, alpha, size): + """ Wrapper around stats.gamma.rvs that converts NegativeBinomial's + parametrization to scipy.gamma. All parameter arrays should have + been broadcasted properly by generate_samples at this point and size is + the scipy.rvs representation. + """ + return stats.gamma.rvs( + a=alpha, + scale=mu / alpha, + size=size, + ) + def logp(self, value): """ Calculate log-probability of ZeroInflatedNegativeBinomial distribution at specified value. diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index ab8988ef9f0..6744f42158a 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -944,3 +944,163 @@ def test_density_dist_without_random_not_sampleable(): samples = 500 with pytest.raises(ValueError): pm.sample_posterior_predictive(trace, samples=samples, model=model, size=100) + + +class TestNestedRandom(SeededTest): + def build_model(self, distribution, shape, nested_rvs_info): + with pm.Model() as model: + nested_rvs = {} + for rv_name, info in nested_rvs_info.items(): + try: + value, nested_shape = info + loc = 0. + except ValueError: + value, nested_shape, loc = info + if value is None: + nested_rvs[rv_name] = pm.Uniform( + rv_name, + 0 + loc, + 1 + loc, + shape=nested_shape, + ) + else: + nested_rvs[rv_name] = value * np.ones(nested_shape) + rv = distribution( + "target", + shape=shape, + **nested_rvs, + ) + return model, rv, nested_rvs + + def sample_prior( + self, + distribution, + shape, + nested_rvs_info, + prior_samples + ): + model, rv, nested_rvs = self.build_model( + distribution, + shape, + nested_rvs_info, + ) + with model: + return pm.sample_prior_predictive(prior_samples) + + @pytest.mark.parametrize( + ["prior_samples", "shape", "psi", "mu", "alpha"], + [ + [10, (3,), (0.5, tuple()), (None, tuple()), (None, (3,))], + [10, (3,), (0.5, (3,)), (None, tuple()), (None, (3,))], + [10, (3,), (0.5, tuple()), (None, (3,)), (None, tuple())], + [10, (3,), (0.5, (3,)), (None, (3,)), (None, tuple())], + [10, (4, 3,), (0.5, (3,)), (None, (3,)), (None, (3,))], + [10, (4, 3,), (0.5, (3,)), (None, (3,)), (None, (4, 3))], + ], + ids=str, + ) + def test_ZeroInflatedNegativeBinomial( + self, + prior_samples, + shape, + psi, + mu, + alpha, + ): + prior = self.sample_prior( + distribution=pm.ZeroInflatedNegativeBinomial, + shape=shape, + nested_rvs_info=dict(psi=psi, mu=mu, alpha=alpha), + prior_samples=prior_samples, + ) + assert prior["target"].shape == (prior_samples,) + shape + + @pytest.mark.parametrize( + ["prior_samples", "shape", "nu", "sigma"], + [ + [10, (3,), (None, tuple()), (None, (3,))], + [10, (3,), (None, tuple()), (None, (3,))], + [10, (3,), (None, (3,)), (None, tuple())], + [10, (3,), (None, (3,)), (None, tuple())], + [10, (4, 3,), (None, (3,)), (None, (3,))], + [10, (4, 3,), (None, (3,)), (None, (4, 3))], + ], + ids=str, + ) + def test_Rice( + self, + prior_samples, + shape, + nu, + sigma, + ): + prior = self.sample_prior( + distribution=pm.Rice, + shape=shape, + nested_rvs_info=dict(nu=nu, sigma=sigma), + prior_samples=prior_samples, + ) + assert prior["target"].shape == (prior_samples,) + shape + + @pytest.mark.parametrize( + ["prior_samples", "shape", "mu", "sigma", "lower", "upper"], + [ + [10, (3,), (None, tuple()), (1., tuple()), (None, tuple()), (None, (3,))], + [10, (3,), (None, tuple()), (1., tuple()), (None, tuple()), (None, (3,))], + [10, (3,), (None, tuple()), (1., tuple()), (None, (3,)), (None, tuple())], + [10, (3,), (None, tuple()), (1., tuple()), (None, (3,)), (None, tuple())], + [10, (4, 3,), (None, (3,)), (1., tuple()), (None, (3,)), (None, (3,))], + [10, (4, 3,), (None, (3,)), (1., tuple()), (None, (3,)), (None, (4, 3))], + [10, (3,), (0., tuple()), (None, tuple()), (None, tuple()), (None, (3,))], + [10, (3,), (0., tuple()), (None, tuple()), (None, tuple()), (None, (3,))], + [10, (3,), (0., tuple()), (None, tuple()), (None, (3,)), (None, tuple())], + [10, (3,), (0., tuple()), (None, tuple()), (None, (3,)), (None, tuple())], + [10, (4, 3,), (0., tuple()), (None, (3,)), (None, (3,)), (None, (3,))], + [10, (4, 3,), (0., tuple()), (None, (3,)), (None, (3,)), (None, (4, 3))], + ], + ids=str, + ) + def test_TruncatedNormal( + self, + prior_samples, + shape, + mu, + sigma, + lower, + upper, + ): + prior = self.sample_prior( + distribution=pm.TruncatedNormal, + shape=shape, + nested_rvs_info=dict(mu=mu, sigma=sigma, lower=lower, upper=upper), + prior_samples=prior_samples, + ) + assert prior["target"].shape == (prior_samples,) + shape + + + @pytest.mark.parametrize( + ["prior_samples", "shape", "c", "lower", "upper"], + [ + [10, (3,), (None, tuple()), (-1., (3,)), (2, tuple())], + [10, (3,), (None, tuple()), (-1., tuple()), (None, tuple(), 1)], + [10, (3,), (None, (3,)), (-1., tuple()), (None, tuple(), 1)], + [10, (4, 3,), (None, (3,)), (-1., tuple()), (None, (3,), 1)], + [10, (4, 3,), (None, (3,)), (None, tuple(), -1), (None, (3,), 1)], + ], + ids=str, + ) + def test_Triangular( + self, + prior_samples, + shape, + c, + lower, + upper, + ): + prior = self.sample_prior( + distribution=pm.Triangular, + shape=shape, + nested_rvs_info=dict(c=c, lower=lower, upper=upper), + prior_samples=prior_samples, + ) + assert prior["target"].shape == (prior_samples,) + shape From 8bfd1a22ba106f44dac8f4657f5e13ad07621c73 Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Thu, 6 Jun 2019 11:24:41 +0200 Subject: [PATCH 2/3] Added release notes --- RELEASE-NOTES.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 3ed800e889a..5447efed802 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -5,6 +5,9 @@ ### New features - Distinguish between `Data` and `Deterministic` variables when graphing models with graphviz. PR [#3491](https://github.com/pymc-defs/pymc3/pulls/3491). +### Maintenance +- Moved math operations out of `Rice`, `TruncatedNormal`, `Triangular` and `ZeroInflatedNegativeBinomial` `random` methods. Math operations on values returned by `draw_values` might not broadcast well, and all the `size` aware broadcasting is left to `generate_samples`. Fixes [#3481](https://github.com/pymc-devs/pymc3/issues/3481) and [#3508](https://github.com/pymc-devs/pymc3/issues/3508) + ## PyMC3 3.7 (May 29 2019) ### New features From da00e555a49d1a4e89df3bea1d79fa1d54fa863d Mon Sep 17 00:00:00 2001 From: lucianopaz Date: Thu, 6 Jun 2019 11:45:22 +0200 Subject: [PATCH 3/3] Fixed domain error in TruncatedNormal test --- pymc3/tests/test_distributions_random.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 6744f42158a..ccfa4937f20 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -1045,18 +1045,18 @@ def test_Rice( @pytest.mark.parametrize( ["prior_samples", "shape", "mu", "sigma", "lower", "upper"], [ - [10, (3,), (None, tuple()), (1., tuple()), (None, tuple()), (None, (3,))], - [10, (3,), (None, tuple()), (1., tuple()), (None, tuple()), (None, (3,))], - [10, (3,), (None, tuple()), (1., tuple()), (None, (3,)), (None, tuple())], - [10, (3,), (None, tuple()), (1., tuple()), (None, (3,)), (None, tuple())], - [10, (4, 3,), (None, (3,)), (1., tuple()), (None, (3,)), (None, (3,))], - [10, (4, 3,), (None, (3,)), (1., tuple()), (None, (3,)), (None, (4, 3))], - [10, (3,), (0., tuple()), (None, tuple()), (None, tuple()), (None, (3,))], - [10, (3,), (0., tuple()), (None, tuple()), (None, tuple()), (None, (3,))], - [10, (3,), (0., tuple()), (None, tuple()), (None, (3,)), (None, tuple())], - [10, (3,), (0., tuple()), (None, tuple()), (None, (3,)), (None, tuple())], - [10, (4, 3,), (0., tuple()), (None, (3,)), (None, (3,)), (None, (3,))], - [10, (4, 3,), (0., tuple()), (None, (3,)), (None, (3,)), (None, (4, 3))], + [10, (3,), (None, tuple()), (1., tuple()), (None, tuple(), -1), (None, (3,))], + [10, (3,), (None, tuple()), (1., tuple()), (None, tuple(), -1), (None, (3,))], + [10, (3,), (None, tuple()), (1., tuple()), (None, (3,), -1), (None, tuple())], + [10, (3,), (None, tuple()), (1., tuple()), (None, (3,), -1), (None, tuple())], + [10, (4, 3,), (None, (3,)), (1., tuple()), (None, (3,), -1), (None, (3,))], + [10, (4, 3,), (None, (3,)), (1., tuple()), (None, (3,), -1), (None, (4, 3))], + [10, (3,), (0., tuple()), (None, tuple()), (None, tuple(), -1), (None, (3,))], + [10, (3,), (0., tuple()), (None, tuple()), (None, tuple(), -1), (None, (3,))], + [10, (3,), (0., tuple()), (None, tuple()), (None, (3,), -1), (None, tuple())], + [10, (3,), (0., tuple()), (None, tuple()), (None, (3,), -1), (None, tuple())], + [10, (4, 3,), (0., tuple()), (None, (3,)), (None, (3,), -1), (None, (3,))], + [10, (4, 3,), (0., tuple()), (None, (3,)), (None, (3,), -1), (None, (4, 3))], ], ids=str, )