From 0e865dc84b3727b0cde0d334fef8cfa226c4dcee Mon Sep 17 00:00:00 2001 From: harivallabha Date: Thu, 17 Sep 2020 22:59:36 +0530 Subject: [PATCH 1/8] Add hypergeometric dist. to discrete --- pymc3/distributions/discrete.py | 107 +++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 1 deletion(-) diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index c5d43706f3..0b505eeec2 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -28,7 +28,8 @@ __all__ = ['Binomial', 'BetaBinomial', 'Bernoulli', 'DiscreteWeibull', 'Poisson', 'NegativeBinomial', 'ConstantDist', 'Constant', 'ZeroInflatedPoisson', 'ZeroInflatedBinomial', 'ZeroInflatedNegativeBinomial', - 'DiscreteUniform', 'Geometric', 'Categorical', 'OrderedLogistic'] + 'DiscreteUniform', 'Geometric', 'HyperGeometric', 'Categorical', + 'OrderedLogistic'] class Binomial(Discrete): @@ -819,6 +820,110 @@ def _repr_latex_(self, name=None, dist=None): get_variable_name(p)) +class HyperGeometric(Discrete): + R""" + Hypergeometric log-likelihood. + + The probability of x successes in a sequence of n Bernoulli + trials (That is, sample size = n) - where the population + size is N, containing a total of k successful individuals. + The process is carried out without replacement. + + The pmf of this distribution is + .. math:: f(x \mid N, n, k) = \frac{\binom{k}{x}\binom{N-k}{n-x}}{\binom{N}{n}} + .. plot:: + import matplotlib.pyplot as plt + import numpy as np + import scipy.stats as st + plt.style.use('seaborn-darkgrid') + x = np.arange(1, 15) + N = 50 + k = 10 + for n in [20, 25]: + pmf = st.hypergeom.pmf(x, N, k, n) + plt.plot(x, pmf, '-o', label='n = {}'.format(n)) + plt.plot(x, pmf, '-o', label='N = {}'.format(N)) + plt.plot(x, pmf, '-o', label='k = {}'.format(k)) + plt.xlabel('x', fontsize=12) + plt.ylabel('f(x)', fontsize=12) + plt.legend(loc=1) + plt.show() + + ======== ============================= + Support :math:`x \in \mathbb{N}_{>0}` + Mean :math:`\dfrac{n.k}{N}` + Variance :math:`\dfrac{(N-n).n.k.(N-k)}{(N-1).N^2}` + ======== ============================= + + Parameters + ---------- + N : integer + Total size of the population + n : integer + Number of samples drawn from the population + k : integer + Number of successful individuals in the population + """ + + def __init__(self, N, k, n, *args, **kwargs): + super().__init__(*args, **kwargs) + self.N = N = tt.as_tensor_variable(intX(N)) + self.k = k = tt.as_tensor_variable(intX(k)) + self.n = n = tt.as_tensor_variable(intX(n)) + self.mode = intX(tt.floor((n + 1)*(k + 1)/(N + 2))) + + def random(self, point=None, size=None): + r""" + Draw random values from HyperGeometric distribution. + Parameters + ---------- + point : dict, optional + Dict of variable values on which random values are to be + conditioned (uses default point if not specified). + size : int, optional + Desired size of random sample (returns one sample if not + specified). + Returns + ------- + array + """ + N, n, k = draw_values([self.N, self.n, self.k], point=point, size=size) + return generate_samples(np.random.hypergeometric, N, n, k, + dist_shape=self.shape, + size=size) + + def logp(self, value): + r""" + Calculate log-probability of HyperGeometric distribution at specified value. + Parameters + ---------- + value : numeric + Value(s) for which log-probability is calculated. If the log probabilities for multiple + values are desired the values must be provided in a numpy array or theano tensor + Returns + ------- + TensorVariable + """ + N = self.N + k = self.k + n = self.n + return bound(binomln(k, value) + binomln(N - k, n - value) - binomln(N, n), + 0 <= k, k <= N, 0 <= n, 0 <= N, n - N + k <= value, 0 <= value, + value <= k, value <= n) + + def _repr_latex_(self, name=None, dist=None): + if dist is None: + dist = self + N = dist.N + n = dist.n + k = dist.k + name = r'\text{%s}' % name + return r'${} \sim \text{{HyperGeometric}}(\mathit{{N}}={},~\mathit{{n}}={},~\mathit{{k}}={})$'.format(name, + get_variable_name(N), + get_variable_name(n), + get_variable_name(k)) + + class DiscreteUniform(Discrete): R""" Discrete uniform distribution. From 7d12a7be96b51782d812aa2803750ce390aceff3 Mon Sep 17 00:00:00 2001 From: harivallabha Date: Thu, 17 Sep 2020 23:00:30 +0530 Subject: [PATCH 2/8] Add tests for the hypergeometric distribution --- pymc3/tests/test_distributions.py | 5 +++++ pymc3/tests/test_distributions_random.py | 7 +++++++ 2 files changed, 12 insertions(+) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index a52ff63a4a..829991d6e1 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -75,6 +75,7 @@ Rice, Kumaraswamy, Moyal, + HyperGeometric ) from ..distributions import continuous @@ -817,6 +818,10 @@ def test_geometric(self): Geometric, Nat, {"p": Unit}, lambda value, p: np.log(sp.geom.pmf(value, p)) ) + def test_hypergeometric(self): + self.pymc3_matches_scipy(HyperGeometric, Nat, {'N': NatSmall, 'n': NatSmall, 'k': NatSmall}, + lambda value, N, n, k: sp.hypergeom.logpmf(value, N, k, n)) + def test_negative_binomial(self): def test_fun(value, mu, alpha): return sp.nbinom.logpmf(value, alpha, 1 - mu / (mu + alpha)) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index b4d4338f1e..1288ce70f0 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -467,6 +467,10 @@ class TestGeometric(BaseTestCases.BaseTestCase): distribution = pm.Geometric params = {'p': 0.5} +class TestHyperGeometric(BaseTestCases.BaseTestCase): + distribution = pm.HyperGeometric + params = {'N': 50, 'n' : 25, 'k' :10} + class TestMoyal(BaseTestCases.BaseTestCase): distribution = pm.Moyal @@ -657,6 +661,9 @@ def ref_rand(size, alpha, mu): def test_geometric(self): pymc3_random_discrete(pm.Geometric, {'p': Unit}, size=500, fails=50, ref_rand=nr.geometric) + def test_hypergeometric(self): + pymc3_random_discrete(pm.HyperGeometric, {'N': Nat, 'n': Nat, 'k': Nat}, size=500, fails=50, ref_rand=nr.hypergeometric) + def test_discrete_uniform(self): def ref_rand(size, lower, upper): return st.randint.rvs(lower, upper + 1, size=size) From a2bffe7e3006582f0d123cecfaac548c3934987f Mon Sep 17 00:00:00 2001 From: harivallabha Date: Thu, 17 Sep 2020 23:05:00 +0530 Subject: [PATCH 3/8] Clean up linting --- pymc3/distributions/__init__.py | 2 ++ pymc3/tests/test_distributions.py | 2 +- pymc3/tests/test_distributions_random.py | 3 ++- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/pymc3/distributions/__init__.py b/pymc3/distributions/__init__.py index bfecad95ef..2b87a59a76 100644 --- a/pymc3/distributions/__init__.py +++ b/pymc3/distributions/__init__.py @@ -63,6 +63,7 @@ from .discrete import ZeroInflatedBinomial from .discrete import DiscreteUniform from .discrete import Geometric +from .discrete import HyperGeometric from .discrete import Categorical from .discrete import OrderedLogistic @@ -136,6 +137,7 @@ 'ZeroInflatedBinomial', 'DiscreteUniform', 'Geometric', + 'HyperGeometric', 'Categorical', 'OrderedLogistic', 'DensityDist', diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 829991d6e1..25d2303c3a 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -821,7 +821,7 @@ def test_geometric(self): def test_hypergeometric(self): self.pymc3_matches_scipy(HyperGeometric, Nat, {'N': NatSmall, 'n': NatSmall, 'k': NatSmall}, lambda value, N, n, k: sp.hypergeom.logpmf(value, N, k, n)) - + def test_negative_binomial(self): def test_fun(value, mu, alpha): return sp.nbinom.logpmf(value, alpha, 1 - mu / (mu + alpha)) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index 1288ce70f0..f52d2065d8 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -662,7 +662,8 @@ def test_geometric(self): pymc3_random_discrete(pm.Geometric, {'p': Unit}, size=500, fails=50, ref_rand=nr.geometric) def test_hypergeometric(self): - pymc3_random_discrete(pm.HyperGeometric, {'N': Nat, 'n': Nat, 'k': Nat}, size=500, fails=50, ref_rand=nr.hypergeometric) + pymc3_random_discrete(pm.HyperGeometric, {'N': Nat, 'n': Nat, 'k': Nat}, size=500, fails=50, + ref_rand=nr.hypergeometric) def test_discrete_uniform(self): def ref_rand(size, lower, upper): From fec3029ddf0a36f62357d8e9972f712221adfbe7 Mon Sep 17 00:00:00 2001 From: harivallabha Date: Thu, 17 Sep 2020 23:08:55 +0530 Subject: [PATCH 4/8] Fix linting - 2 --- pymc3/tests/test_distributions_random.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index f52d2065d8..c0b54c1162 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -467,6 +467,7 @@ class TestGeometric(BaseTestCases.BaseTestCase): distribution = pm.Geometric params = {'p': 0.5} + class TestHyperGeometric(BaseTestCases.BaseTestCase): distribution = pm.HyperGeometric params = {'N': 50, 'n' : 25, 'k' :10} From a23960fb93c91cb6ad21888b458150941192c754 Mon Sep 17 00:00:00 2001 From: harivallabha Date: Thu, 17 Sep 2020 23:31:10 +0530 Subject: [PATCH 5/8] Add line in RELEASE-NOTES.md --- RELEASE-NOTES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index ac313f2877..ed24de4128 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -11,6 +11,7 @@ ### New features - `sample_posterior_predictive_w` can now feed on `xarray.Dataset` - e.g. from `InferenceData.posterior`. (see [#4042](https://github.com/pymc-devs/pymc3/pull/4042)) +- Support HyperGeometric Distribution through `pymc3.distributions.discrete.HyperGeometric`. (see [#4108](https://github.com/pymc-devs/pymc3/pull/4108)) ## PyMC3 3.9.3 (11 August 2020) From bd6e7fb876b2f0bd0e355eb12a8f451fb14187b2 Mon Sep 17 00:00:00 2001 From: harivallabha Date: Fri, 18 Sep 2020 01:19:21 +0530 Subject: [PATCH 6/8] Run black formatter on discrete.py, test_distributions.py and test_distributions_random.py --- pymc3/distributions/discrete.py | 423 ++++---- pymc3/tests/test_distributions.py | 19 +- pymc3/tests/test_distributions_random.py | 1133 ++++++++++++++-------- 3 files changed, 988 insertions(+), 587 deletions(-) diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index 0b505eeec2..919d75c059 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -25,11 +25,24 @@ from ..theanof import floatX, intX, take_along_axis -__all__ = ['Binomial', 'BetaBinomial', 'Bernoulli', 'DiscreteWeibull', - 'Poisson', 'NegativeBinomial', 'ConstantDist', 'Constant', - 'ZeroInflatedPoisson', 'ZeroInflatedBinomial', 'ZeroInflatedNegativeBinomial', - 'DiscreteUniform', 'Geometric', 'HyperGeometric', 'Categorical', - 'OrderedLogistic'] +__all__ = [ + "Binomial", + "BetaBinomial", + "Bernoulli", + "DiscreteWeibull", + "Poisson", + "NegativeBinomial", + "ConstantDist", + "Constant", + "ZeroInflatedPoisson", + "ZeroInflatedBinomial", + "ZeroInflatedNegativeBinomial", + "DiscreteUniform", + "Geometric", + "HyperGeometric", + "Categorical", + "OrderedLogistic", +] class Binomial(Discrete): @@ -98,9 +111,9 @@ def random(self, point=None, size=None): array """ n, p = draw_values([self.n, self.p], point=point, size=size) - return generate_samples(stats.binom.rvs, n=n, p=p, - dist_shape=self.shape, - size=size) + return generate_samples( + stats.binom.rvs, n=n, p=p, dist_shape=self.shape, size=size + ) def logp(self, value): r""" @@ -121,18 +134,22 @@ def logp(self, value): return bound( binomln(n, value) + logpow(p, value) + logpow(1 - p, n - value), - 0 <= value, value <= n, - 0 <= p, p <= 1) + 0 <= value, + value <= n, + 0 <= p, + p <= 1, + ) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self n = dist.n p = dist.p - name = r'\text{%s}' % name - return r'${} \sim \text{{Binomial}}(\mathit{{n}}={},~\mathit{{p}}={})$'.format(name, - get_variable_name(n), - get_variable_name(p)) + name = r"\text{%s}" % name + return r"${} \sim \text{{Binomial}}(\mathit{{n}}={},~\mathit{{p}}={})$".format( + name, get_variable_name(n), get_variable_name(p) + ) + class BetaBinomial(Discrete): R""" @@ -194,7 +211,7 @@ def __init__(self, alpha, beta, n, *args, **kwargs): self.alpha = alpha = tt.as_tensor_variable(floatX(alpha)) self.beta = beta = tt.as_tensor_variable(floatX(beta)) self.n = n = tt.as_tensor_variable(intX(n)) - self.mode = tt.cast(tround(alpha / (alpha + beta)), 'int8') + self.mode = tt.cast(tround(alpha / (alpha + beta)), "int8") def _random(self, alpha, beta, n, size=None): size = size or 1 @@ -208,8 +225,11 @@ def _random(self, alpha, beta, n, size=None): quotient, remainder = divmod(_p.shape[0], _n.shape[0]) if remainder != 0: - raise TypeError('n has a bad size! Was cast to {}, must evenly divide {}'.format( - _n.shape[0], _p.shape[0])) + raise TypeError( + "n has a bad size! Was cast to {}, must evenly divide {}".format( + _n.shape[0], _p.shape[0] + ) + ) if quotient != 1: _n = np.tile(_n, quotient) samples = np.reshape(stats.binom.rvs(n=_n, p=_p, size=_size), size) @@ -232,11 +252,12 @@ def random(self, point=None, size=None): ------- array """ - alpha, beta, n = \ - draw_values([self.alpha, self.beta, self.n], point=point, size=size) - return generate_samples(self._random, alpha=alpha, beta=beta, n=n, - dist_shape=self.shape, - size=size) + alpha, beta, n = draw_values( + [self.alpha, self.beta, self.n], point=point, size=size + ) + return generate_samples( + self._random, alpha=alpha, beta=beta, n=n, dist_shape=self.shape, size=size + ) def logp(self, value): r""" @@ -254,21 +275,25 @@ def logp(self, value): """ alpha = self.alpha beta = self.beta - return bound(binomln(self.n, value) - + betaln(value + alpha, self.n - value + beta) - - betaln(alpha, beta), - value >= 0, value <= self.n, - alpha > 0, beta > 0) + return bound( + binomln(self.n, value) + + betaln(value + alpha, self.n - value + beta) + - betaln(alpha, beta), + value >= 0, + value <= self.n, + alpha > 0, + beta > 0, + ) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self alpha = dist.alpha beta = dist.beta - name = r'\text{%s}' % name - return r'${} \sim \text{{BetaBinomial}}(\mathit{{alpha}}={},~\mathit{{beta}}={})$'.format(name, - get_variable_name(alpha), - get_variable_name(beta)) + name = r"\text{%s}" % name + return r"${} \sim \text{{BetaBinomial}}(\mathit{{alpha}}={},~\mathit{{beta}}={})$".format( + name, get_variable_name(alpha), get_variable_name(beta) + ) class Bernoulli(Discrete): @@ -314,7 +339,7 @@ class Bernoulli(Discrete): def __init__(self, p=None, logit_p=None, *args, **kwargs): super().__init__(*args, **kwargs) if sum(int(var is None) for var in [p, logit_p]) != 1: - raise ValueError('Specify one of p and logit_p') + raise ValueError("Specify one of p and logit_p") if p is not None: self._is_logit = False self.p = p = tt.as_tensor_variable(floatX(p)) @@ -324,7 +349,7 @@ def __init__(self, p=None, logit_p=None, *args, **kwargs): self.p = tt.nnet.sigmoid(floatX(logit_p)) self._logit_p = tt.as_tensor_variable(logit_p) - self.mode = tt.cast(tround(self.p), 'int8') + self.mode = tt.cast(tround(self.p), "int8") def random(self, point=None, size=None): r""" @@ -344,9 +369,9 @@ def random(self, point=None, size=None): array """ p = draw_values([self.p], point=point, size=size)[0] - return generate_samples(stats.bernoulli.rvs, p, - dist_shape=self.shape, - size=size) + return generate_samples( + stats.bernoulli.rvs, p, dist_shape=self.shape, size=size + ) def logp(self, value): r""" @@ -369,16 +394,20 @@ def logp(self, value): p = self.p return bound( tt.switch(value, tt.log(p), tt.log(1 - p)), - value >= 0, value <= 1, - p >= 0, p <= 1) + value >= 0, + value <= 1, + p >= 0, + p <= 1, + ) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self p = dist.p - name = r'\text{%s}' % name - return r'${} \sim \text{{Bernoulli}}(\mathit{{p}}={})$'.format(name, - get_variable_name(p)) + name = r"\text{%s}" % name + return r"${} \sim \text{{Bernoulli}}(\mathit{{p}}={})$".format( + name, get_variable_name(p) + ) class DiscreteWeibull(Discrete): @@ -419,8 +448,9 @@ def DiscreteWeibull(q, b, x): Variance :math:`2 \sum_{x = 1}^{\infty} x q^{x^{\beta}} - \mu - \mu^2` ======== ====================== """ + def __init__(self, q, beta, *args, **kwargs): - super().__init__(*args, defaults=('median',), **kwargs) + super().__init__(*args, defaults=("median",), **kwargs) self.q = q = tt.as_tensor_variable(floatX(q)) self.beta = beta = tt.as_tensor_variable(floatX(beta)) @@ -444,10 +474,16 @@ def logp(self, value): q = self.q beta = self.beta - return bound(tt.log(tt.power(q, tt.power(value, beta)) - tt.power(q, tt.power(value + 1, beta))), - 0 <= value, - 0 < q, q < 1, - 0 < beta) + return bound( + tt.log( + tt.power(q, tt.power(value, beta)) + - tt.power(q, tt.power(value + 1, beta)) + ), + 0 <= value, + 0 < q, + q < 1, + 0 < beta, + ) def _ppf(self, p): r""" @@ -457,12 +493,14 @@ def _ppf(self, p): q = self.q beta = self.beta - return (tt.ceil(tt.power(tt.log(1 - p) / tt.log(q), 1. / beta)) - 1).astype('int64') + return (tt.ceil(tt.power(tt.log(1 - p) / tt.log(q), 1.0 / beta)) - 1).astype( + "int64" + ) def _random(self, q, beta, size=None): p = np.random.uniform(size=size) - return np.ceil(np.power(np.log(1 - p) / np.log(q), 1. / beta)) - 1 + return np.ceil(np.power(np.log(1 - p) / np.log(q), 1.0 / beta)) - 1 def random(self, point=None, size=None): r""" @@ -483,19 +521,17 @@ def random(self, point=None, size=None): """ q, beta = draw_values([self.q, self.beta], point=point, size=size) - return generate_samples(self._random, q, beta, - dist_shape=self.shape, - size=size) + return generate_samples(self._random, q, beta, dist_shape=self.shape, size=size) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self q = dist.q beta = dist.beta - name = r'\text{%s}' % name - return r'${} \sim \text{{DiscreteWeibull}}(\mathit{{q}}={},~\mathit{{beta}}={})$'.format(name, - get_variable_name(q), - get_variable_name(beta)) + name = r"\text{%s}" % name + return r"${} \sim \text{{DiscreteWeibull}}(\mathit{{q}}={},~\mathit{{beta}}={})$".format( + name, get_variable_name(q), get_variable_name(beta) + ) class Poisson(Discrete): @@ -565,9 +601,7 @@ def random(self, point=None, size=None): array """ mu = draw_values([self.mu], point=point, size=size)[0] - return generate_samples(stats.poisson.rvs, mu, - dist_shape=self.shape, - size=size) + return generate_samples(stats.poisson.rvs, mu, dist_shape=self.shape, size=size) def logp(self, value): r""" @@ -584,20 +618,18 @@ def logp(self, value): TensorVariable """ mu = self.mu - log_prob = bound( - logpow(mu, value) - factln(value) - mu, - mu >= 0, value >= 0) + log_prob = bound(logpow(mu, value) - factln(value) - mu, mu >= 0, value >= 0) # Return zero when mu and value are both zero - return tt.switch(tt.eq(mu, 0) * tt.eq(value, 0), - 0, log_prob) + return tt.switch(tt.eq(mu, 0) * tt.eq(value, 0), 0, log_prob) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self mu = dist.mu - name = r'\text{%s}' % name - return r'${} \sim \text{{Poisson}}(\mathit{{mu}}={})$'.format(name, - get_variable_name(mu)) + name = r"\text{%s}" % name + return r"${} \sim \text{{Poisson}}(\mathit{{mu}}={})$".format( + name, get_variable_name(mu) + ) class NegativeBinomial(Discrete): @@ -674,14 +706,14 @@ def random(self, point=None, size=None): array """ mu, alpha = draw_values([self.mu, self.alpha], point=point, size=size) - g = generate_samples(self._random, mu=mu, alpha=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 return np.asarray(stats.poisson.rvs(g)).reshape(g.shape) def _random(self, mu, alpha, size): - r""" Wrapper around stats.gamma.rvs that converts NegativeBinomial's + r"""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. @@ -708,25 +740,29 @@ def logp(self, value): """ mu = self.mu alpha = self.alpha - negbinom = bound(binomln(value + alpha - 1, value) - + logpow(mu / (mu + alpha), value) - + logpow(alpha / (mu + alpha), alpha), - value >= 0, mu > 0, alpha > 0) + negbinom = bound( + binomln(value + alpha - 1, value) + + logpow(mu / (mu + alpha), value) + + logpow(alpha / (mu + alpha), alpha), + value >= 0, + mu > 0, + alpha > 0, + ) # Return Poisson when alpha gets very large. - return tt.switch(tt.gt(alpha, 1e10), - Poisson.dist(self.mu).logp(value), - negbinom) + return tt.switch( + tt.gt(alpha, 1e10), Poisson.dist(self.mu).logp(value), negbinom + ) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self mu = dist.mu alpha = dist.alpha - name = r'\text{%s}' % name - return r'${} \sim \text{{NegativeBinomial}}(\mathit{{mu}}={},~\mathit{{alpha}}={})$'.format(name, - get_variable_name(mu), - get_variable_name(alpha)) + name = r"\text{%s}" % name + return r"${} \sim \text{{NegativeBinomial}}(\mathit{{mu}}={},~\mathit{{alpha}}={})$".format( + name, get_variable_name(mu), get_variable_name(alpha) + ) class Geometric(Discrete): @@ -789,9 +825,9 @@ def random(self, point=None, size=None): array """ p = draw_values([self.p], point=point, size=size)[0] - return generate_samples(np.random.geometric, p, - dist_shape=self.shape, - size=size) + return generate_samples( + np.random.geometric, p, dist_shape=self.shape, size=size + ) def logp(self, value): r""" @@ -808,16 +844,16 @@ def logp(self, value): TensorVariable """ p = self.p - return bound(tt.log(p) + logpow(1 - p, value - 1), - 0 <= p, p <= 1, value >= 1) + return bound(tt.log(p) + logpow(1 - p, value - 1), 0 <= p, p <= 1, value >= 1) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self p = dist.p - name = r'\text{%s}' % name - return r'${} \sim \text{{Geometric}}(\mathit{{p}}={})$'.format(name, - get_variable_name(p)) + name = r"\text{%s}" % name + return r"${} \sim \text{{Geometric}}(\mathit{{p}}={})$".format( + name, get_variable_name(p) + ) class HyperGeometric(Discrete): @@ -865,12 +901,12 @@ class HyperGeometric(Discrete): Number of successful individuals in the population """ - def __init__(self, N, k, n, *args, **kwargs): + def __init__(self, N, k, n, *args, **kwargs): super().__init__(*args, **kwargs) self.N = N = tt.as_tensor_variable(intX(N)) self.k = k = tt.as_tensor_variable(intX(k)) self.n = n = tt.as_tensor_variable(intX(n)) - self.mode = intX(tt.floor((n + 1)*(k + 1)/(N + 2))) + self.mode = intX(tt.floor((n + 1) * (k + 1) / (N + 2))) def random(self, point=None, size=None): r""" @@ -888,9 +924,9 @@ def random(self, point=None, size=None): array """ N, n, k = draw_values([self.N, self.n, self.k], point=point, size=size) - return generate_samples(np.random.hypergeometric, N, n, k, - dist_shape=self.shape, - size=size) + return generate_samples( + np.random.hypergeometric, N, n, k, dist_shape=self.shape, size=size + ) def logp(self, value): r""" @@ -907,9 +943,17 @@ def logp(self, value): N = self.N k = self.k n = self.n - return bound(binomln(k, value) + binomln(N - k, n - value) - binomln(N, n), - 0 <= k, k <= N, 0 <= n, 0 <= N, n - N + k <= value, 0 <= value, - value <= k, value <= n) + return bound( + binomln(k, value) + binomln(N - k, n - value) - binomln(N, n), + 0 <= k, + k <= N, + 0 <= n, + 0 <= N, + n - N + k <= value, + 0 <= value, + value <= k, + value <= n, + ) def _repr_latex_(self, name=None, dist=None): if dist is None: @@ -917,11 +961,10 @@ def _repr_latex_(self, name=None, dist=None): N = dist.N n = dist.n k = dist.k - name = r'\text{%s}' % name - return r'${} \sim \text{{HyperGeometric}}(\mathit{{N}}={},~\mathit{{n}}={},~\mathit{{k}}={})$'.format(name, - get_variable_name(N), - get_variable_name(n), - get_variable_name(k)) + name = r"\text{%s}" % name + return r"${} \sim \text{{HyperGeometric}}(\mathit{{N}}={},~\mathit{{n}}={},~\mathit{{k}}={})$".format( + name, get_variable_name(N), get_variable_name(n), get_variable_name(k) + ) class DiscreteUniform(Discrete): @@ -967,8 +1010,7 @@ def __init__(self, lower, upper, *args, **kwargs): super().__init__(*args, **kwargs) self.lower = intX(tt.floor(lower)) self.upper = intX(tt.floor(upper)) - self.mode = tt.maximum( - intX(tt.floor((upper + lower) / 2.)), self.lower) + self.mode = tt.maximum(intX(tt.floor((upper + lower) / 2.0)), self.lower) def _random(self, lower, upper, size=None): # This way seems to be the only to deal with lower and upper @@ -994,10 +1036,9 @@ def random(self, point=None, size=None): array """ lower, upper = draw_values([self.lower, self.upper], point=point, size=size) - return generate_samples(self._random, - lower, upper, - dist_shape=self.shape, - size=size) + return generate_samples( + self._random, lower, upper, dist_shape=self.shape, size=size + ) def logp(self, value): r""" @@ -1015,18 +1056,17 @@ def logp(self, value): """ upper = self.upper lower = self.lower - return bound(-tt.log(upper - lower + 1), - lower <= value, value <= upper) + return bound(-tt.log(upper - lower + 1), lower <= value, value <= upper) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self lower = dist.lower upper = dist.upper - name = r'\text{%s}' % name - return r'${} \sim \text{{DiscreteUniform}}(\mathit{{lower}}={},~\mathit{{upper}}={})$'.format(name, - get_variable_name(lower), - get_variable_name(upper)) + name = r"\text{%s}" % name + return r"${} \sim \text{{DiscreteUniform}}(\mathit{{lower}}={},~\mathit{{upper}}={})$".format( + name, get_variable_name(lower), get_variable_name(upper) + ) class Categorical(Discrete): @@ -1099,11 +1139,13 @@ def random(self, point=None, size=None): p, k = draw_values([self.p, self.k], point=point, size=size) p = p / np.sum(p, axis=-1, keepdims=True) - return generate_samples(random_choice, - p=p, - broadcast_shape=p.shape[:-1] or (1,), - dist_shape=self.shape, - size=size) + return generate_samples( + random_choice, + p=p, + broadcast_shape=p.shape[:-1] or (1,), + dist_shape=self.shape, + size=size, + ) def logp(self, value): r""" @@ -1129,13 +1171,9 @@ def logp(self, value): if p.ndim > 1: if p.ndim > value_clip.ndim: - value_clip = tt.shape_padleft( - value_clip, p_.ndim - value_clip.ndim - ) + value_clip = tt.shape_padleft(value_clip, p_.ndim - value_clip.ndim) elif p.ndim < value_clip.ndim: - p = tt.shape_padleft( - p, value_clip.ndim - p_.ndim - ) + p = tt.shape_padleft(p, value_clip.ndim - p_.ndim) pattern = (p.ndim - 1,) + tuple(range(p.ndim - 1)) a = tt.log( take_along_axis( @@ -1146,16 +1184,22 @@ def logp(self, value): else: a = tt.log(p[value_clip]) - return bound(a, value >= 0, value <= (k - 1), - tt.all(p_ >= 0, axis=-1), tt.all(p <= 1, axis=-1)) + return bound( + a, + value >= 0, + value <= (k - 1), + tt.all(p_ >= 0, axis=-1), + tt.all(p <= 1, axis=-1), + ) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self p = dist.p - name = r'\text{%s}' % name - return r'${} \sim \text{{Categorical}}(\mathit{{p}}={})$'.format(name, - get_variable_name(p)) + name = r"\text{%s}" % name + return r"${} \sim \text{{Categorical}}(\mathit{{p}}={})$".format( + name, get_variable_name(p) + ) class Constant(Discrete): @@ -1169,8 +1213,10 @@ class Constant(Discrete): """ def __init__(self, c, *args, **kwargs): - warnings.warn("Constant has been deprecated. We recommend using a Deterministic object instead.", - DeprecationWarning) + warnings.warn( + "Constant has been deprecated. We recommend using a Deterministic object instead.", + DeprecationWarning, + ) super().__init__(*args, **kwargs) self.mean = self.median = self.mode = self.c = c = tt.as_tensor_variable(c) @@ -1197,8 +1243,9 @@ def random(self, point=None, size=None): def _random(c, dtype=dtype, size=None): return np.full(size, fill_value=c, dtype=dtype) - return generate_samples(_random, c=c, dist_shape=self.shape, - size=size).astype(dtype) + return generate_samples(_random, c=c, dist_shape=self.shape, size=size).astype( + dtype + ) def logp(self, value): r""" @@ -1220,8 +1267,8 @@ def logp(self, value): def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self - name = r'\text{%s}' % name - return r'${} \sim \text{{Constant}}()$'.format(name) + name = r"\text{%s}" % name + return r"${} \sim \text{{Constant}}()$".format(name) ConstantDist = Constant @@ -1302,9 +1349,7 @@ def random(self, point=None, size=None): array """ theta, psi = draw_values([self.theta, self.psi], point=point, size=size) - g = generate_samples(stats.poisson.rvs, theta, - dist_shape=self.shape, - size=size) + g = generate_samples(stats.poisson.rvs, theta, dist_shape=self.shape, size=size) g, psi = broadcast_distribution_samples([g, psi], size=size) return g * (np.random.random(g.shape) < psi) @@ -1328,23 +1373,20 @@ def logp(self, value): logp_val = tt.switch( tt.gt(value, 0), tt.log(psi) + self.pois.logp(value), - logaddexp(tt.log1p(-psi), tt.log(psi) - theta)) + logaddexp(tt.log1p(-psi), tt.log(psi) - theta), + ) - return bound( - logp_val, - 0 <= value, - 0 <= psi, psi <= 1, - 0 <= theta) + return bound(logp_val, 0 <= value, 0 <= psi, psi <= 1, 0 <= theta) def _repr_latex_(self, name=None, dist=None): if dist is None: dist = self theta = dist.theta psi = dist.psi - name = r'\text{%s}' % name - return r'${} \sim \text{{ZeroInflatedPoisson}}(\mathit{{theta}}={},~\mathit{{psi}}={})$'.format(name, - get_variable_name(theta), - get_variable_name(psi)) + name = r"\text{%s}" % name + return r"${} \sim \text{{ZeroInflatedPoisson}}(\mathit{{theta}}={},~\mathit{{psi}}={})$".format( + name, get_variable_name(theta), get_variable_name(psi) + ) class ZeroInflatedBinomial(Discrete): @@ -1424,9 +1466,7 @@ def random(self, point=None, size=None): array """ n, p, psi = draw_values([self.n, self.p, self.psi], point=point, size=size) - g = generate_samples(stats.binom.rvs, n, p, - dist_shape=self.shape, - size=size) + g = generate_samples(stats.binom.rvs, n, p, dist_shape=self.shape, size=size) g, psi = broadcast_distribution_samples([g, psi], size=size) return g * (np.random.random(g.shape) < psi) @@ -1451,13 +1491,12 @@ def logp(self, value): logp_val = tt.switch( tt.gt(value, 0), tt.log(psi) + self.bin.logp(value), - logaddexp(tt.log1p(-psi), tt.log(psi) + n * tt.log1p(-p))) + logaddexp(tt.log1p(-psi), tt.log(psi) + n * tt.log1p(-p)), + ) return bound( - logp_val, - 0 <= value, value <= n, - 0 <= psi, psi <= 1, - 0 <= p, p <= 1) + logp_val, 0 <= value, value <= n, 0 <= psi, psi <= 1, 0 <= p, p <= 1 + ) def _repr_latex_(self, name=None, dist=None): if dist is None: @@ -1469,11 +1508,12 @@ def _repr_latex_(self, name=None, dist=None): name_n = get_variable_name(n) name_p = get_variable_name(p) name_psi = get_variable_name(psi) - name = r'\text{%s}' % name - return (r'${} \sim \text{{ZeroInflatedBinomial}}' - r'(\mathit{{n}}={},~\mathit{{p}}={},~' - r'\mathit{{psi}}={})$' - .format(name, name_n, name_p, name_psi)) + name = r"\text{%s}" % name + return ( + r"${} \sim \text{{ZeroInflatedBinomial}}" + r"(\mathit{{n}}={},~\mathit{{p}}={},~" + r"\mathit{{psi}}={})$".format(name, name_n, name_p, name_psi) + ) class ZeroInflatedNegativeBinomial(Discrete): @@ -1570,20 +1610,17 @@ def random(self, point=None, size=None): array """ mu, alpha, psi = draw_values( - [self.mu, self.alpha, self.psi], point=point, size=size) + [self.mu, self.alpha, self.psi], point=point, size=size + ) g = generate_samples( - self._random, - mu=mu, - alpha=alpha, - dist_shape=self.shape, - size=size + 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): - r""" Wrapper around stats.gamma.rvs that converts NegativeBinomial's + r"""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. @@ -1614,19 +1651,12 @@ def logp(self, value): logp_other = tt.log(psi) + self.nb.logp(value) logp_0 = logaddexp( - tt.log1p(-psi), - tt.log(psi) + alpha * (tt.log(alpha) - tt.log(alpha + mu))) + tt.log1p(-psi), tt.log(psi) + alpha * (tt.log(alpha) - tt.log(alpha + mu)) + ) - logp_val = tt.switch( - tt.gt(value, 0), - logp_other, - logp_0) + logp_val = tt.switch(tt.gt(value, 0), logp_other, logp_0) - return bound( - logp_val, - 0 <= value, - 0 <= psi, psi <= 1, - mu > 0, alpha > 0) + return bound(logp_val, 0 <= value, 0 <= psi, psi <= 1, mu > 0, alpha > 0) def _repr_latex_(self, name=None, dist=None): if dist is None: @@ -1638,11 +1668,12 @@ def _repr_latex_(self, name=None, dist=None): name_mu = get_variable_name(mu) name_alpha = get_variable_name(alpha) name_psi = get_variable_name(psi) - name = r'\text{%s}' % name - return (r'${} \sim \text{{ZeroInflatedNegativeBinomial}}' - r'(\mathit{{mu}}={},~\mathit{{alpha}}={},~' - r'\mathit{{psi}}={})$' - .format(name, name_mu, name_alpha, name_psi)) + name = r"\text{%s}" % name + return ( + r"${} \sim \text{{ZeroInflatedNegativeBinomial}}" + r"(\mathit{{mu}}={},~\mathit{{alpha}}={},~" + r"\mathit{{psi}}={})$".format(name, name_mu, name_alpha, name_psi) + ) class OrderedLogistic(Categorical): @@ -1716,11 +1747,14 @@ def __init__(self, eta, cutpoints, *args, **kwargs): self.cutpoints = tt.as_tensor_variable(cutpoints) pa = sigmoid(self.cutpoints - tt.shape_padright(self.eta)) - p_cum = tt.concatenate([ - tt.zeros_like(tt.shape_padright(pa[..., 0])), - pa, - tt.ones_like(tt.shape_padright(pa[..., 0])) - ], axis=-1) + p_cum = tt.concatenate( + [ + tt.zeros_like(tt.shape_padright(pa[..., 0])), + pa, + tt.ones_like(tt.shape_padright(pa[..., 0])), + ], + axis=-1, + ) p = p_cum[..., 1:] - p_cum[..., :-1] super().__init__(p=p, *args, **kwargs) @@ -1730,6 +1764,9 @@ def _repr_latex_(self, name=None, dist=None): dist = self name_eta = get_variable_name(dist.eta) name_cutpoints = get_variable_name(dist.cutpoints) - return (r'${} \sim \text{{OrderedLogistic}}' - r'(\mathit{{eta}}={}, \mathit{{cutpoints}}={}$' - .format(name, name_eta, name_cutpoints)) + return ( + r"${} \sim \text{{OrderedLogistic}}" + r"(\mathit{{eta}}={}, \mathit{{cutpoints}}={}$".format( + name, name_eta, name_cutpoints + ) + ) diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index 25d2303c3a..5d42794f61 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -75,7 +75,7 @@ Rice, Kumaraswamy, Moyal, - HyperGeometric + HyperGeometric, ) from ..distributions import continuous @@ -819,8 +819,12 @@ def test_geometric(self): ) def test_hypergeometric(self): - self.pymc3_matches_scipy(HyperGeometric, Nat, {'N': NatSmall, 'n': NatSmall, 'k': NatSmall}, - lambda value, N, n, k: sp.hypergeom.logpmf(value, N, k, n)) + self.pymc3_matches_scipy( + HyperGeometric, + Nat, + {"N": NatSmall, "n": NatSmall, "k": NatSmall}, + lambda value, N, n, k: sp.hypergeom.logpmf(value, N, k, n), + ) def test_negative_binomial(self): def test_fun(value, mu, alpha): @@ -1340,7 +1344,9 @@ def test_dirichlet_shape(self): dir_rv = Dirichlet.dist(a) assert dir_rv.shape == (2,) - with pytest.warns(DeprecationWarning), theano.change_flags(compute_test_value="ignore"): + with pytest.warns(DeprecationWarning), theano.change_flags( + compute_test_value="ignore" + ): dir_rv = Dirichlet.dist(tt.vector()) def test_dirichlet_2D(self): @@ -1885,9 +1891,10 @@ def func(x): return -2 * (x ** 2).sum() with pm.Model(): - pm.Normal('x') - y = pm.DensityDist('y', func) + pm.Normal("x") + y = pm.DensityDist("y", func) pm.sample(draws=5, tune=1, mp_ctx="spawn") import pickle + pickle.loads(pickle.dumps(y)) diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py index c0b54c1162..f3d8148bf3 100644 --- a/pymc3/tests/test_distributions_random.py +++ b/pymc3/tests/test_distributions_random.py @@ -23,20 +23,46 @@ import pymc3 as pm from pymc3.distributions.dist_math import clipped_beta_rvs -from pymc3.distributions.distribution import (draw_values, - _DrawValuesContext, - _DrawValuesContextBlocker) +from pymc3.distributions.distribution import ( + draw_values, + _DrawValuesContext, + _DrawValuesContextBlocker, +) from .helpers import SeededTest from .test_distributions import ( - build_model, Domain, product, R, Rplus, Rplusbig, Runif, Rplusdunif, - Unit, Nat, NatSmall, I, Simplex, Vector, PdMatrix, - PdMatrixChol, PdMatrixCholUpper, RealMatrix, RandomPdMatrix + build_model, + Domain, + product, + R, + Rplus, + Rplusbig, + Runif, + Rplusdunif, + Unit, + Nat, + NatSmall, + I, + Simplex, + Vector, + PdMatrix, + PdMatrixChol, + PdMatrixCholUpper, + RealMatrix, + RandomPdMatrix, ) -def pymc3_random(dist, paramdomains, ref_rand, valuedomain=Domain([0]), - size=10000, alpha=0.05, fails=10, extra_args=None, - model_args=None): +def pymc3_random( + dist, + paramdomains, + ref_rand, + valuedomain=Domain([0]), + size=10000, + alpha=0.05, + fails=10, + extra_args=None, + model_args=None, +): if model_args is None: model_args = {} model = build_model(dist, valuedomain, paramdomains, extra_args) @@ -49,17 +75,22 @@ def pymc3_random(dist, paramdomains, ref_rand, valuedomain=Domain([0]), # a certain number of times. Crude, but necessary. f = fails while p <= alpha and f > 0: - s0 = model.named_vars['value'].random(size=size, point=pt) + s0 = model.named_vars["value"].random(size=size, point=pt) s1 = ref_rand(size=size, **pt) - _, p = st.ks_2samp(np.atleast_1d(s0).flatten(), - np.atleast_1d(s1).flatten()) + _, p = st.ks_2samp(np.atleast_1d(s0).flatten(), np.atleast_1d(s1).flatten()) f -= 1 assert p > alpha, str(pt) -def pymc3_random_discrete(dist, paramdomains, - valuedomain=Domain([0]), ref_rand=None, - size=100000, alpha=0.05, fails=20): +def pymc3_random_discrete( + dist, + paramdomains, + valuedomain=Domain([0]), + ref_rand=None, + size=100000, + alpha=0.05, + fails=20, +): model = build_model(dist, valuedomain, paramdomains) domains = paramdomains.copy() for pt in product(domains, n_samples=100): @@ -69,7 +100,7 @@ def pymc3_random_discrete(dist, paramdomains, # a certain number of times. f = fails while p <= alpha and f > 0: - o = model.named_vars['value'].random(size=size, point=pt) + o = model.named_vars["value"].random(size=size, point=pt) e = ref_rand(size=size, **pt) o = np.atleast_1d(o).flatten() e = np.atleast_1d(e).flatten() @@ -79,7 +110,7 @@ def pymc3_random_discrete(dist, paramdomains, expected[e] = (observed.get(e, 0), expected[e]) k = np.array([v for v in expected.values()]) if np.all(k[:, 0] == k[:, 1]): - p = 1. + p = 1.0 else: _, p = st.chisquare(k[:, 0], k[:, 1]) f -= 1 @@ -89,23 +120,23 @@ def pymc3_random_discrete(dist, paramdomains, class TestDrawValues(SeededTest): def test_draw_scalar_parameters(self): with pm.Model(): - y = pm.Normal('y1', mu=0., sigma=1.) + y = pm.Normal("y1", mu=0.0, sigma=1.0) mu, tau = draw_values([y.distribution.mu, y.distribution.tau]) npt.assert_almost_equal(mu, 0) npt.assert_almost_equal(tau, 1) def test_draw_dependencies(self): with pm.Model(): - x = pm.Normal('x', mu=0., sigma=1.) - exp_x = pm.Deterministic('exp_x', pm.math.exp(x)) + x = pm.Normal("x", mu=0.0, sigma=1.0) + exp_x = pm.Deterministic("exp_x", pm.math.exp(x)) x, exp_x = draw_values([x, exp_x]) npt.assert_almost_equal(np.exp(x), exp_x) def test_draw_order(self): with pm.Model(): - x = pm.Normal('x', mu=0., sigma=1.) - exp_x = pm.Deterministic('exp_x', pm.math.exp(x)) + x = pm.Normal("x", mu=0.0, sigma=1.0) + exp_x = pm.Deterministic("exp_x", pm.math.exp(x)) # Need to draw x before drawing log_x exp_x, x = draw_values([exp_x, x]) @@ -113,19 +144,20 @@ def test_draw_order(self): def test_draw_point_replacement(self): with pm.Model(): - mu = pm.Normal('mu', mu=0., tau=1e-3) - sigma = pm.Gamma('sigma', alpha=1., beta=1., transform=None) - y = pm.Normal('y', mu=mu, sigma=sigma) - mu2, tau2 = draw_values([y.distribution.mu, y.distribution.tau], - point={'mu': 5., 'sigma': 2.}) + mu = pm.Normal("mu", mu=0.0, tau=1e-3) + sigma = pm.Gamma("sigma", alpha=1.0, beta=1.0, transform=None) + y = pm.Normal("y", mu=mu, sigma=sigma) + mu2, tau2 = draw_values( + [y.distribution.mu, y.distribution.tau], point={"mu": 5.0, "sigma": 2.0} + ) npt.assert_almost_equal(mu2, 5) - npt.assert_almost_equal(tau2, 1 / 2.**2) + npt.assert_almost_equal(tau2, 1 / 2.0 ** 2) def test_random_sample_returns_nd_array(self): with pm.Model(): - mu = pm.Normal('mu', mu=0., tau=1e-3) - sigma = pm.Gamma('sigma', alpha=1., beta=1., transform=None) - y = pm.Normal('y', mu=mu, sigma=sigma) + mu = pm.Normal("mu", mu=0.0, tau=1e-3) + sigma = pm.Gamma("sigma", alpha=1.0, beta=1.0, transform=None) + y = pm.Normal("y", mu=mu, sigma=sigma) mu, tau = draw_values([y.distribution.mu, y.distribution.tau]) assert isinstance(mu, np.ndarray) assert isinstance(tau, np.ndarray) @@ -135,38 +167,38 @@ class TestDrawValuesContext: def test_normal_context(self): with _DrawValuesContext() as context0: assert context0.parent is None - context0.drawn_vars['root_test'] = 1 + context0.drawn_vars["root_test"] = 1 with _DrawValuesContext() as context1: assert id(context1.drawn_vars) == id(context0.drawn_vars) assert context1.parent == context0 with _DrawValuesContext() as context2: assert id(context2.drawn_vars) == id(context0.drawn_vars) assert context2.parent == context1 - context2.drawn_vars['leaf_test'] = 2 - assert context1.drawn_vars['leaf_test'] == 2 - context1.drawn_vars['root_test'] = 3 - assert context0.drawn_vars['root_test'] == 3 - assert context0.drawn_vars['leaf_test'] == 2 + context2.drawn_vars["leaf_test"] = 2 + assert context1.drawn_vars["leaf_test"] == 2 + context1.drawn_vars["root_test"] = 3 + assert context0.drawn_vars["root_test"] == 3 + assert context0.drawn_vars["leaf_test"] == 2 def test_blocking_context(self): with _DrawValuesContext() as context0: assert context0.parent is None - context0.drawn_vars['root_test'] = 1 + context0.drawn_vars["root_test"] = 1 with _DrawValuesContext() as context1: assert id(context1.drawn_vars) == id(context0.drawn_vars) assert context1.parent == context0 with _DrawValuesContextBlocker() as blocker: assert id(blocker.drawn_vars) != id(context0.drawn_vars) assert blocker.parent is None - blocker.drawn_vars['root_test'] = 2 + blocker.drawn_vars["root_test"] = 2 with _DrawValuesContext() as context2: assert id(context2.drawn_vars) == id(blocker.drawn_vars) assert context2.parent == blocker - context2.drawn_vars['root_test'] = 3 - context2.drawn_vars['leaf_test'] = 4 - assert blocker.drawn_vars['root_test'] == 3 - assert 'leaf_test' not in context1.drawn_vars - assert context0.drawn_vars['root_test'] == 1 + context2.drawn_vars["root_test"] = 3 + context2.drawn_vars["leaf_test"] = 4 + assert blocker.drawn_vars["root_test"] == 3 + assert "leaf_test" not in context1.drawn_vars + assert context0.drawn_vars["root_test"] == 1 class BaseTestCases: @@ -179,8 +211,10 @@ def setup_method(self, *args, **kwargs): def get_random_variable(self, shape, with_vector_params=False, name=None): if with_vector_params: - params = {key: value * np.ones(self.shape, dtype=np.dtype(type(value))) for - key, value in self.params.items()} + params = { + key: value * np.ones(self.shape, dtype=np.dtype(type(value))) + for key, value in self.params.items() + } else: params = self.params if name is None: @@ -190,7 +224,9 @@ def get_random_variable(self, shape, with_vector_params=False, name=None): return self.distribution(name, transform=None, **params) else: try: - return self.distribution(name, shape=shape, transform=None, **params) + return self.distribution( + name, shape=shape, transform=None, **params + ) except TypeError: if np.sum(np.atleast_1d(shape)) == 0: pytest.skip("Timeseries must have positive shape") @@ -203,17 +239,17 @@ def sample_random_variable(random_variable, size): except AttributeError: return random_variable.distribution.random(size=size) - @pytest.mark.parametrize('size', [None, 5, (4, 5)], ids=str) + @pytest.mark.parametrize("size", [None, 5, (4, 5)], ids=str) def test_scalar_parameter_shape(self, size): rv = self.get_random_variable(None) if size is None: - expected = 1, + expected = (1,) else: expected = np.atleast_1d(size).tolist() actual = np.atleast_1d(self.sample_random_variable(rv, size)).shape assert tuple(expected) == actual - @pytest.mark.parametrize('size', [None, 5, (4, 5)], ids=str) + @pytest.mark.parametrize("size", [None, 5, (4, 5)], ids=str) def test_scalar_shape(self, size): shape = 10 rv = self.get_random_variable(shape) @@ -226,7 +262,7 @@ def test_scalar_shape(self, size): actual = np.atleast_1d(self.sample_random_variable(rv, size)).shape assert tuple(expected) == actual - @pytest.mark.parametrize('size', [None, 5, (4, 5)], ids=str) + @pytest.mark.parametrize("size", [None, 5, (4, 5)], ids=str) def test_parameters_1d_shape(self, size): rv = self.get_random_variable(self.shape, with_vector_params=True) if size is None: @@ -237,7 +273,7 @@ def test_parameters_1d_shape(self, size): actual = self.sample_random_variable(rv, size).shape assert tuple(expected) == actual - @pytest.mark.parametrize('size', [None, 5, (4, 5)], ids=str) + @pytest.mark.parametrize("size", [None, 5, (4, 5)], ids=str) def test_broadcast_shape(self, size): broadcast_shape = (2 * self.shape, self.shape) rv = self.get_random_variable(broadcast_shape, with_vector_params=True) @@ -249,11 +285,13 @@ def test_broadcast_shape(self, size): actual = np.atleast_1d(self.sample_random_variable(rv, size)).shape assert tuple(expected) == actual - @pytest.mark.parametrize('shape', [(), (1,), (1, 1), (1, 2), (10, 10, 1), (10, 10, 2)], ids=str) + @pytest.mark.parametrize( + "shape", [(), (1,), (1, 1), (1, 2), (10, 10, 1), (10, 10, 2)], ids=str + ) def test_different_shapes_and_sample_sizes(self, shape): prefix = self.distribution.__name__ - rv = self.get_random_variable(shape, name='%s_%s' % (prefix, shape)) + rv = self.get_random_variable(shape, name="%s_%s" % (prefix, shape)) for size in (None, 1, 5, (4, 5)): if size is None: s = [] @@ -273,216 +311,225 @@ def test_different_shapes_and_sample_sizes(self, shape): class TestGaussianRandomWalk(BaseTestCases.BaseTestCase): distribution = pm.GaussianRandomWalk - params = {'mu': 1., 'sigma': 1.} + params = {"mu": 1.0, "sigma": 1.0} @pytest.mark.xfail(reason="Supporting this makes a nasty API") def test_broadcast_shape(self): super().test_broadcast_shape() + class TestNormal(BaseTestCases.BaseTestCase): distribution = pm.Normal - params = {'mu': 0., 'tau': 1.} + params = {"mu": 0.0, "tau": 1.0} + class TestTruncatedNormal(BaseTestCases.BaseTestCase): distribution = pm.TruncatedNormal - params = {'mu': 0., 'tau': 1., 'lower': -0.5, 'upper': 0.5} + params = {"mu": 0.0, "tau": 1.0, "lower": -0.5, "upper": 0.5} + class TestTruncatedNormalLower(BaseTestCases.BaseTestCase): distribution = pm.TruncatedNormal - params = {'mu': 0., 'tau': 1., 'lower': -0.5} + params = {"mu": 0.0, "tau": 1.0, "lower": -0.5} + class TestTruncatedNormalUpper(BaseTestCases.BaseTestCase): distribution = pm.TruncatedNormal - params = {'mu': 0., 'tau': 1., 'upper': 0.5} + params = {"mu": 0.0, "tau": 1.0, "upper": 0.5} + class TestSkewNormal(BaseTestCases.BaseTestCase): distribution = pm.SkewNormal - params = {'mu': 0., 'sigma': 1., 'alpha': 5.} + params = {"mu": 0.0, "sigma": 1.0, "alpha": 5.0} class TestHalfNormal(BaseTestCases.BaseTestCase): distribution = pm.HalfNormal - params = {'tau': 1.} + params = {"tau": 1.0} class TestUniform(BaseTestCases.BaseTestCase): distribution = pm.Uniform - params = {'lower': 0., 'upper': 1.} + params = {"lower": 0.0, "upper": 1.0} class TestTriangular(BaseTestCases.BaseTestCase): distribution = pm.Triangular - params = {'c': 0.5, 'lower': 0., 'upper': 1.} + params = {"c": 0.5, "lower": 0.0, "upper": 1.0} class TestWald(BaseTestCases.BaseTestCase): distribution = pm.Wald - params = {'mu': 1., 'lam': 1., 'alpha': 0.} + params = {"mu": 1.0, "lam": 1.0, "alpha": 0.0} class TestBeta(BaseTestCases.BaseTestCase): distribution = pm.Beta - params = {'alpha': 1., 'beta': 1.} + params = {"alpha": 1.0, "beta": 1.0} class TestKumaraswamy(BaseTestCases.BaseTestCase): distribution = pm.Kumaraswamy - params = {'a': 1., 'b': 1.} + params = {"a": 1.0, "b": 1.0} class TestExponential(BaseTestCases.BaseTestCase): distribution = pm.Exponential - params = {'lam': 1.} + params = {"lam": 1.0} class TestLaplace(BaseTestCases.BaseTestCase): distribution = pm.Laplace - params = {'mu': 1., 'b': 1.} + params = {"mu": 1.0, "b": 1.0} class TestLognormal(BaseTestCases.BaseTestCase): distribution = pm.Lognormal - params = {'mu': 1., 'tau': 1.} + params = {"mu": 1.0, "tau": 1.0} class TestStudentT(BaseTestCases.BaseTestCase): distribution = pm.StudentT - params = {'nu': 5., 'mu': 0., 'lam': 1.} + params = {"nu": 5.0, "mu": 0.0, "lam": 1.0} class TestPareto(BaseTestCases.BaseTestCase): distribution = pm.Pareto - params = {'alpha': 0.5, 'm': 1.} + params = {"alpha": 0.5, "m": 1.0} class TestCauchy(BaseTestCases.BaseTestCase): distribution = pm.Cauchy - params = {'alpha': 1., 'beta': 1.} + params = {"alpha": 1.0, "beta": 1.0} class TestHalfCauchy(BaseTestCases.BaseTestCase): distribution = pm.HalfCauchy - params = {'beta': 1.} + params = {"beta": 1.0} class TestGamma(BaseTestCases.BaseTestCase): distribution = pm.Gamma - params = {'alpha': 1., 'beta': 1.} + params = {"alpha": 1.0, "beta": 1.0} class TestInverseGamma(BaseTestCases.BaseTestCase): distribution = pm.InverseGamma - params = {'alpha': 0.5, 'beta': 0.5} + params = {"alpha": 0.5, "beta": 0.5} class TestChiSquared(BaseTestCases.BaseTestCase): distribution = pm.ChiSquared - params = {'nu': 2.} + params = {"nu": 2.0} class TestWeibull(BaseTestCases.BaseTestCase): distribution = pm.Weibull - params = {'alpha': 1., 'beta': 1.} + params = {"alpha": 1.0, "beta": 1.0} class TestExGaussian(BaseTestCases.BaseTestCase): distribution = pm.ExGaussian - params = {'mu': 0., 'sigma': 1., 'nu': 1.} + params = {"mu": 0.0, "sigma": 1.0, "nu": 1.0} class TestVonMises(BaseTestCases.BaseTestCase): distribution = pm.VonMises - params = {'mu': 0., 'kappa': 1.} + params = {"mu": 0.0, "kappa": 1.0} class TestGumbel(BaseTestCases.BaseTestCase): distribution = pm.Gumbel - params = {'mu': 0., 'beta': 1.} + params = {"mu": 0.0, "beta": 1.0} class TestLogistic(BaseTestCases.BaseTestCase): distribution = pm.Logistic - params = {'mu': 0., 's': 1.} + params = {"mu": 0.0, "s": 1.0} class TestLogitNormal(BaseTestCases.BaseTestCase): distribution = pm.LogitNormal - params = {'mu': 0., 'sigma': 1.} + params = {"mu": 0.0, "sigma": 1.0} class TestBinomial(BaseTestCases.BaseTestCase): distribution = pm.Binomial - params = {'n': 5, 'p': 0.5} + params = {"n": 5, "p": 0.5} class TestBetaBinomial(BaseTestCases.BaseTestCase): distribution = pm.BetaBinomial - params = {'n': 5, 'alpha': 1., 'beta': 1.} + params = {"n": 5, "alpha": 1.0, "beta": 1.0} class TestBernoulli(BaseTestCases.BaseTestCase): distribution = pm.Bernoulli - params = {'p': 0.5} + params = {"p": 0.5} class TestDiscreteWeibull(BaseTestCases.BaseTestCase): distribution = pm.DiscreteWeibull - params = {'q': 0.25, 'beta': 2.} + params = {"q": 0.25, "beta": 2.0} class TestPoisson(BaseTestCases.BaseTestCase): distribution = pm.Poisson - params = {'mu': 1.} + params = {"mu": 1.0} class TestNegativeBinomial(BaseTestCases.BaseTestCase): distribution = pm.NegativeBinomial - params = {'mu': 1., 'alpha': 1.} + params = {"mu": 1.0, "alpha": 1.0} class TestConstant(BaseTestCases.BaseTestCase): distribution = pm.Constant - params = {'c': 3} + params = {"c": 3} class TestZeroInflatedPoisson(BaseTestCases.BaseTestCase): distribution = pm.ZeroInflatedPoisson - params = {'theta': 1., 'psi': 0.3} + params = {"theta": 1.0, "psi": 0.3} class TestZeroInflatedNegativeBinomial(BaseTestCases.BaseTestCase): distribution = pm.ZeroInflatedNegativeBinomial - params = {'mu': 1., 'alpha': 1., 'psi': 0.3} + params = {"mu": 1.0, "alpha": 1.0, "psi": 0.3} + class TestZeroInflatedBinomial(BaseTestCases.BaseTestCase): distribution = pm.ZeroInflatedBinomial - params = {'n': 10, 'p': 0.6, 'psi': 0.3} + params = {"n": 10, "p": 0.6, "psi": 0.3} + class TestDiscreteUniform(BaseTestCases.BaseTestCase): distribution = pm.DiscreteUniform - params = {'lower': 0., 'upper': 10.} + params = {"lower": 0.0, "upper": 10.0} class TestGeometric(BaseTestCases.BaseTestCase): distribution = pm.Geometric - params = {'p': 0.5} + params = {"p": 0.5} class TestHyperGeometric(BaseTestCases.BaseTestCase): distribution = pm.HyperGeometric - params = {'N': 50, 'n' : 25, 'k' :10} + params = {"N": 50, "n": 25, "k": 10} + - class TestMoyal(BaseTestCases.BaseTestCase): distribution = pm.Moyal - params = {'mu': 0., 'sigma': 1.} + params = {"mu": 0.0, "sigma": 1.0} + - class TestCategorical(BaseTestCases.BaseTestCase): distribution = pm.Categorical - params = {'p': np.ones(BaseTestCases.BaseTestCase.shape)} + params = {"p": np.ones(BaseTestCases.BaseTestCase.shape)} - def get_random_variable(self, shape, with_vector_params=False, **kwargs): # don't transform categories + def get_random_variable( + self, shape, with_vector_params=False, **kwargs + ): # don't transform categories return super().get_random_variable(shape, with_vector_params=False, **kwargs) def test_probability_vector_shape(self): @@ -502,195 +549,284 @@ def test_bounded(self): def ref_rand(size, tau): return -st.halfnorm.rvs(size=size, loc=0, scale=tau ** -0.5) - pymc3_random(BoundedNormal, {'tau': Rplus}, ref_rand=ref_rand) + + pymc3_random(BoundedNormal, {"tau": Rplus}, ref_rand=ref_rand) def test_uniform(self): def ref_rand(size, lower, upper): return st.uniform.rvs(size=size, loc=lower, scale=upper - lower) - pymc3_random(pm.Uniform, {'lower': -Rplus, 'upper': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.Uniform, {"lower": -Rplus, "upper": Rplus}, ref_rand=ref_rand) def test_normal(self): def ref_rand(size, mu, sigma): return st.norm.rvs(size=size, loc=mu, scale=sigma) - pymc3_random(pm.Normal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.Normal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand) def test_truncated_normal(self): def ref_rand(size, mu, sigma, lower, upper): - return st.truncnorm.rvs((lower - mu) / sigma, (upper - mu) / sigma, size=size, loc=mu, scale=sigma) - pymc3_random(pm.TruncatedNormal, {'mu': R, 'sigma': Rplusbig, 'lower': -Rplusbig, 'upper': Rplusbig}, - ref_rand=ref_rand) + return st.truncnorm.rvs( + (lower - mu) / sigma, + (upper - mu) / sigma, + size=size, + loc=mu, + scale=sigma, + ) + + pymc3_random( + pm.TruncatedNormal, + {"mu": R, "sigma": Rplusbig, "lower": -Rplusbig, "upper": Rplusbig}, + ref_rand=ref_rand, + ) def test_truncated_normal_lower(self): def ref_rand(size, mu, sigma, lower): - return st.truncnorm.rvs((lower - mu) / sigma, np.inf, size=size, loc=mu, scale=sigma) - pymc3_random(pm.TruncatedNormal, {'mu': R, 'sigma': Rplusbig, 'lower': -Rplusbig}, - ref_rand=ref_rand) + return st.truncnorm.rvs( + (lower - mu) / sigma, np.inf, size=size, loc=mu, scale=sigma + ) + + pymc3_random( + pm.TruncatedNormal, + {"mu": R, "sigma": Rplusbig, "lower": -Rplusbig}, + ref_rand=ref_rand, + ) def test_truncated_normal_upper(self): def ref_rand(size, mu, sigma, upper): - return st.truncnorm.rvs(-np.inf, (upper - mu) / sigma, size=size, loc=mu, scale=sigma) - pymc3_random(pm.TruncatedNormal, {'mu': R, 'sigma': Rplusbig, 'upper': Rplusbig}, - ref_rand=ref_rand) + return st.truncnorm.rvs( + -np.inf, (upper - mu) / sigma, size=size, loc=mu, scale=sigma + ) + + pymc3_random( + pm.TruncatedNormal, + {"mu": R, "sigma": Rplusbig, "upper": Rplusbig}, + ref_rand=ref_rand, + ) def test_skew_normal(self): def ref_rand(size, alpha, mu, sigma): return st.skewnorm.rvs(size=size, a=alpha, loc=mu, scale=sigma) - pymc3_random(pm.SkewNormal, {'mu': R, 'sigma': Rplus, 'alpha': R}, ref_rand=ref_rand) + + pymc3_random( + pm.SkewNormal, {"mu": R, "sigma": Rplus, "alpha": R}, ref_rand=ref_rand + ) def test_half_normal(self): def ref_rand(size, tau): return st.halfnorm.rvs(size=size, loc=0, scale=tau ** -0.5) - pymc3_random(pm.HalfNormal, {'tau': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.HalfNormal, {"tau": Rplus}, ref_rand=ref_rand) def test_wald(self): # Cannot do anything too exciting as scipy wald is a # location-scale model of the *standard* wald with mu=1 and lam=1 def ref_rand(size, mu, lam, alpha): return st.wald.rvs(size=size, loc=alpha) - pymc3_random(pm.Wald, - {'mu': Domain([1., 1., 1.]), 'lam': Domain( - [1., 1., 1.]), 'alpha': Rplus}, - ref_rand=ref_rand) + + pymc3_random( + pm.Wald, + { + "mu": Domain([1.0, 1.0, 1.0]), + "lam": Domain([1.0, 1.0, 1.0]), + "alpha": Rplus, + }, + ref_rand=ref_rand, + ) def test_beta(self): def ref_rand(size, alpha, beta): return clipped_beta_rvs(a=alpha, b=beta, size=size) - pymc3_random(pm.Beta, {'alpha': Rplus, 'beta': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.Beta, {"alpha": Rplus, "beta": Rplus}, ref_rand=ref_rand) def test_exponential(self): def ref_rand(size, lam): - return nr.exponential(scale=1. / lam, size=size) - pymc3_random(pm.Exponential, {'lam': Rplus}, ref_rand=ref_rand) + return nr.exponential(scale=1.0 / lam, size=size) + + pymc3_random(pm.Exponential, {"lam": Rplus}, ref_rand=ref_rand) def test_laplace(self): def ref_rand(size, mu, b): return st.laplace.rvs(mu, b, size=size) - pymc3_random(pm.Laplace, {'mu': R, 'b': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.Laplace, {"mu": R, "b": Rplus}, ref_rand=ref_rand) def test_lognormal(self): def ref_rand(size, mu, tau): - return np.exp(mu + (tau ** -0.5) * st.norm.rvs(loc=0., scale=1., size=size)) - pymc3_random(pm.Lognormal, {'mu': R, 'tau': Rplusbig}, ref_rand=ref_rand) + return np.exp( + mu + (tau ** -0.5) * st.norm.rvs(loc=0.0, scale=1.0, size=size) + ) + + pymc3_random(pm.Lognormal, {"mu": R, "tau": Rplusbig}, ref_rand=ref_rand) def test_student_t(self): def ref_rand(size, nu, mu, lam): - return st.t.rvs(nu, mu, lam**-.5, size=size) - pymc3_random(pm.StudentT, {'nu': Rplus, 'mu': R, 'lam': Rplus}, ref_rand=ref_rand) + return st.t.rvs(nu, mu, lam ** -0.5, size=size) + + pymc3_random( + pm.StudentT, {"nu": Rplus, "mu": R, "lam": Rplus}, ref_rand=ref_rand + ) def test_cauchy(self): def ref_rand(size, alpha, beta): return st.cauchy.rvs(alpha, beta, size=size) - pymc3_random(pm.Cauchy, {'alpha': R, 'beta': Rplusbig}, ref_rand=ref_rand) + + pymc3_random(pm.Cauchy, {"alpha": R, "beta": Rplusbig}, ref_rand=ref_rand) def test_half_cauchy(self): def ref_rand(size, beta): return st.halfcauchy.rvs(scale=beta, size=size) - pymc3_random(pm.HalfCauchy, {'beta': Rplusbig}, ref_rand=ref_rand) + + pymc3_random(pm.HalfCauchy, {"beta": Rplusbig}, ref_rand=ref_rand) def test_gamma_alpha_beta(self): def ref_rand(size, alpha, beta): - return st.gamma.rvs(alpha, scale=1. / beta, size=size) - pymc3_random(pm.Gamma, {'alpha': Rplusbig, 'beta': Rplusbig}, ref_rand=ref_rand) + return st.gamma.rvs(alpha, scale=1.0 / beta, size=size) + + pymc3_random(pm.Gamma, {"alpha": Rplusbig, "beta": Rplusbig}, ref_rand=ref_rand) def test_gamma_mu_sigma(self): def ref_rand(size, mu, sigma): - return st.gamma.rvs(mu**2 / sigma**2, scale=sigma ** 2 / mu, size=size) - pymc3_random(pm.Gamma, {'mu': Rplusbig, 'sigma': Rplusbig}, ref_rand=ref_rand) + return st.gamma.rvs(mu ** 2 / sigma ** 2, scale=sigma ** 2 / mu, size=size) + + pymc3_random(pm.Gamma, {"mu": Rplusbig, "sigma": Rplusbig}, ref_rand=ref_rand) def test_inverse_gamma(self): def ref_rand(size, alpha, beta): return st.invgamma.rvs(a=alpha, scale=beta, size=size) - pymc3_random(pm.InverseGamma, {'alpha': Rplus, 'beta': Rplus}, ref_rand=ref_rand) + + pymc3_random( + pm.InverseGamma, {"alpha": Rplus, "beta": Rplus}, ref_rand=ref_rand + ) def test_pareto(self): def ref_rand(size, alpha, m): return st.pareto.rvs(alpha, scale=m, size=size) - pymc3_random(pm.Pareto, {'alpha': Rplusbig, 'm': Rplusbig}, ref_rand=ref_rand) + + pymc3_random(pm.Pareto, {"alpha": Rplusbig, "m": Rplusbig}, ref_rand=ref_rand) def test_ex_gaussian(self): def ref_rand(size, mu, sigma, nu): return nr.normal(mu, sigma, size=size) + nr.exponential(scale=nu, size=size) - pymc3_random(pm.ExGaussian, {'mu': R, 'sigma': Rplus, 'nu': Rplus}, ref_rand=ref_rand) + + pymc3_random( + pm.ExGaussian, {"mu": R, "sigma": Rplus, "nu": Rplus}, ref_rand=ref_rand + ) def test_vonmises(self): def ref_rand(size, mu, kappa): return st.vonmises.rvs(size=size, loc=mu, kappa=kappa) - pymc3_random(pm.VonMises, {'mu': R, 'kappa': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.VonMises, {"mu": R, "kappa": Rplus}, ref_rand=ref_rand) def test_triangular(self): def ref_rand(size, lower, upper, c): scale = upper - lower c_ = (c - lower) / scale return st.triang.rvs(size=size, loc=lower, scale=scale, c=c_) - pymc3_random(pm.Triangular, {'lower': Runif, 'upper': Runif + 3, 'c': Runif + 1}, ref_rand=ref_rand) + + pymc3_random( + pm.Triangular, + {"lower": Runif, "upper": Runif + 3, "c": Runif + 1}, + ref_rand=ref_rand, + ) def test_flat(self): with pm.Model(): - f = pm.Flat('f') + f = pm.Flat("f") with pytest.raises(ValueError): f.random(1) def test_half_flat(self): with pm.Model(): - f = pm.HalfFlat('f') + f = pm.HalfFlat("f") with pytest.raises(ValueError): f.random(1) def test_binomial(self): - pymc3_random_discrete(pm.Binomial, {'n': Nat, 'p': Unit}, ref_rand=st.binom.rvs) + pymc3_random_discrete(pm.Binomial, {"n": Nat, "p": Unit}, ref_rand=st.binom.rvs) def test_beta_binomial(self): - pymc3_random_discrete(pm.BetaBinomial, {'n': Nat, 'alpha': Rplus, 'beta': Rplus}, - ref_rand=self._beta_bin) + pymc3_random_discrete( + pm.BetaBinomial, + {"n": Nat, "alpha": Rplus, "beta": Rplus}, + ref_rand=self._beta_bin, + ) def _beta_bin(self, n, alpha, beta, size=None): return st.binom.rvs(n, st.beta.rvs(a=alpha, b=beta, size=size)) def test_bernoulli(self): - pymc3_random_discrete(pm.Bernoulli, {'p': Unit}, - ref_rand=lambda size, p=None: st.bernoulli.rvs(p, size=size)) + pymc3_random_discrete( + pm.Bernoulli, + {"p": Unit}, + ref_rand=lambda size, p=None: st.bernoulli.rvs(p, size=size), + ) def test_poisson(self): - pymc3_random_discrete(pm.Poisson, {'mu': Rplusbig}, size=500, ref_rand=st.poisson.rvs) + pymc3_random_discrete( + pm.Poisson, {"mu": Rplusbig}, size=500, ref_rand=st.poisson.rvs + ) def test_negative_binomial(self): def ref_rand(size, alpha, mu): return st.nbinom.rvs(alpha, alpha / (mu + alpha), size=size) - pymc3_random_discrete(pm.NegativeBinomial, {'mu': Rplusbig, 'alpha': Rplusbig}, - size=100, fails=50, ref_rand=ref_rand) + + pymc3_random_discrete( + pm.NegativeBinomial, + {"mu": Rplusbig, "alpha": Rplusbig}, + size=100, + fails=50, + ref_rand=ref_rand, + ) def test_geometric(self): - pymc3_random_discrete(pm.Geometric, {'p': Unit}, size=500, fails=50, ref_rand=nr.geometric) + pymc3_random_discrete( + pm.Geometric, {"p": Unit}, size=500, fails=50, ref_rand=nr.geometric + ) def test_hypergeometric(self): - pymc3_random_discrete(pm.HyperGeometric, {'N': Nat, 'n': Nat, 'k': Nat}, size=500, fails=50, - ref_rand=nr.hypergeometric) + pymc3_random_discrete( + pm.HyperGeometric, + {"N": Nat, "n": Nat, "k": Nat}, + size=500, + fails=50, + ref_rand=nr.hypergeometric, + ) def test_discrete_uniform(self): def ref_rand(size, lower, upper): return st.randint.rvs(lower, upper + 1, size=size) - pymc3_random_discrete(pm.DiscreteUniform, {'lower': -NatSmall, 'upper': NatSmall}, - ref_rand=ref_rand) + + pymc3_random_discrete( + pm.DiscreteUniform, + {"lower": -NatSmall, "upper": NatSmall}, + ref_rand=ref_rand, + ) def test_discrete_weibull(self): def ref_rand(size, q, beta): u = np.random.uniform(size=size) - return np.ceil(np.power(np.log(1 - u) / np.log(q), 1. / beta)) - 1 + return np.ceil(np.power(np.log(1 - u) / np.log(q), 1.0 / beta)) - 1 - pymc3_random_discrete(pm.DiscreteWeibull, {'q': Unit, 'beta': Rplusdunif}, - ref_rand=ref_rand) + pymc3_random_discrete( + pm.DiscreteWeibull, {"q": Unit, "beta": Rplusdunif}, ref_rand=ref_rand + ) - @pytest.mark.parametrize('s', [2, 3, 4]) + @pytest.mark.parametrize("s", [2, 3, 4]) def test_categorical_random(self, s): def ref_rand(size, p): return nr.choice(np.arange(p.shape[0]), p=p, size=size) - pymc3_random_discrete(pm.Categorical, {'p': Simplex(s)}, ref_rand=ref_rand) + + pymc3_random_discrete(pm.Categorical, {"p": Simplex(s)}, ref_rand=ref_rand) def test_constant_dist(self): def ref_rand(size, c): return c * np.ones(size, dtype=int) - pymc3_random_discrete(pm.Constant, {'c': I}, ref_rand=ref_rand) + + pymc3_random_discrete(pm.Constant, {"c": I}, ref_rand=ref_rand) def test_mv_normal(self): def ref_rand(size, mu, cov): @@ -706,41 +842,82 @@ def ref_rand_uchol(size, mu, chol): return ref_rand(size, mu, np.dot(chol.T, chol)) for n in [2, 3]: - pymc3_random(pm.MvNormal, {'mu': Vector(R, n), 'cov': PdMatrix(n)}, - size=100, valuedomain=Vector(R, n), ref_rand=ref_rand) - pymc3_random(pm.MvNormal, {'mu': Vector(R, n), 'tau': PdMatrix(n)}, - size=100, valuedomain=Vector(R, n), ref_rand=ref_rand_tau) - pymc3_random(pm.MvNormal, {'mu': Vector(R, n), 'chol': PdMatrixChol(n)}, - size=100, valuedomain=Vector(R, n), ref_rand=ref_rand_chol) pymc3_random( pm.MvNormal, - {'mu': Vector(R, n), 'chol': PdMatrixCholUpper(n)}, - size=100, valuedomain=Vector(R, n), ref_rand=ref_rand_uchol, - extra_args={'lower': False} + {"mu": Vector(R, n), "cov": PdMatrix(n)}, + size=100, + valuedomain=Vector(R, n), + ref_rand=ref_rand, + ) + pymc3_random( + pm.MvNormal, + {"mu": Vector(R, n), "tau": PdMatrix(n)}, + size=100, + valuedomain=Vector(R, n), + ref_rand=ref_rand_tau, + ) + pymc3_random( + pm.MvNormal, + {"mu": Vector(R, n), "chol": PdMatrixChol(n)}, + size=100, + valuedomain=Vector(R, n), + ref_rand=ref_rand_chol, + ) + pymc3_random( + pm.MvNormal, + {"mu": Vector(R, n), "chol": PdMatrixCholUpper(n)}, + size=100, + valuedomain=Vector(R, n), + ref_rand=ref_rand_uchol, + extra_args={"lower": False}, ) def test_matrix_normal(self): def ref_rand(size, mu, rowcov, colcov): - return st.matrix_normal.rvs(mean=mu, rowcov=rowcov, colcov=colcov, size=size) + return st.matrix_normal.rvs( + mean=mu, rowcov=rowcov, colcov=colcov, size=size + ) # def ref_rand_tau(size, mu, tau): # return ref_rand(size, mu, linalg.inv(tau)) def ref_rand_chol(size, mu, rowchol, colchol): - return ref_rand(size, mu, rowcov=np.dot(rowchol, rowchol.T), - colcov=np.dot(colchol, colchol.T)) + return ref_rand( + size, + mu, + rowcov=np.dot(rowchol, rowchol.T), + colcov=np.dot(colchol, colchol.T), + ) def ref_rand_uchol(size, mu, rowchol, colchol): - return ref_rand(size, mu, rowcov=np.dot(rowchol.T, rowchol), - colcov=np.dot(colchol.T, colchol)) + return ref_rand( + size, + mu, + rowcov=np.dot(rowchol.T, rowchol), + colcov=np.dot(colchol.T, colchol), + ) for n in [2, 3]: - pymc3_random(pm.MatrixNormal, {'mu': RealMatrix(n, n), 'rowcov': PdMatrix(n), 'colcov': PdMatrix(n)}, - size=n, valuedomain=RealMatrix(n, n), ref_rand=ref_rand) + pymc3_random( + pm.MatrixNormal, + {"mu": RealMatrix(n, n), "rowcov": PdMatrix(n), "colcov": PdMatrix(n)}, + size=n, + valuedomain=RealMatrix(n, n), + ref_rand=ref_rand, + ) # pymc3_random(pm.MatrixNormal, {'mu': RealMatrix(n, n), 'tau': PdMatrix(n)}, # size=n, valuedomain=RealMatrix(n, n), ref_rand=ref_rand_tau) - pymc3_random(pm.MatrixNormal, {'mu': RealMatrix(n, n), 'rowchol': PdMatrixChol(n), 'colchol': PdMatrixChol(n)}, - size=n, valuedomain=RealMatrix(n, n), ref_rand=ref_rand_chol) + pymc3_random( + pm.MatrixNormal, + { + "mu": RealMatrix(n, n), + "rowchol": PdMatrixChol(n), + "colchol": PdMatrixChol(n), + }, + size=n, + valuedomain=RealMatrix(n, n), + ref_rand=ref_rand_chol, + ) # pymc3_random( # pm.MvNormal, # {'mu': RealMatrix(n, n), 'rowchol': PdMatrixCholUpper(n), 'colchol': PdMatrixCholUpper(n)}, @@ -751,7 +928,7 @@ def ref_rand_uchol(size, mu, rowchol, colchol): def test_kronecker_normal(self): def ref_rand(size, mu, covs, sigma): cov = pm.math.kronecker(covs[0], covs[1]).eval() - cov += sigma**2 * np.identity(cov.shape[0]) + cov += sigma ** 2 * np.identity(cov.shape[0]) return st.multivariate_normal.rvs(mean=mu, cov=cov, size=size) def ref_rand_chol(size, mu, chols, sigma): @@ -767,104 +944,141 @@ def ref_rand_evd(size, mu, evds, sigma): sizes = [2, 3] sigmas = [0, 1] for n, sigma in zip(sizes, sigmas): - N = n**2 + N = n ** 2 covs = [RandomPdMatrix(n), RandomPdMatrix(n)] chols = list(map(np.linalg.cholesky, covs)) evds = list(map(np.linalg.eigh, covs)) - dom = Domain([np.random.randn(N)*0.1], edges=(None, None), shape=N) - mu = Domain([np.random.randn(N)*0.1], edges=(None, None), shape=N) + dom = Domain([np.random.randn(N) * 0.1], edges=(None, None), shape=N) + mu = Domain([np.random.randn(N) * 0.1], edges=(None, None), shape=N) - std_args = {'mu': mu} - cov_args = {'covs': covs} - chol_args = {'chols': chols} - evd_args = {'evds': evds} + std_args = {"mu": mu} + cov_args = {"covs": covs} + chol_args = {"chols": chols} + evd_args = {"evds": evds} if sigma is not None and sigma != 0: - std_args['sigma'] = Domain([sigma], edges=(None, None)) + std_args["sigma"] = Domain([sigma], edges=(None, None)) else: for args in [cov_args, chol_args, evd_args]: - args['sigma'] = sigma + args["sigma"] = sigma pymc3_random( - pm.KroneckerNormal, std_args, valuedomain=dom, - ref_rand=ref_rand, extra_args=cov_args, model_args=cov_args) + pm.KroneckerNormal, + std_args, + valuedomain=dom, + ref_rand=ref_rand, + extra_args=cov_args, + model_args=cov_args, + ) pymc3_random( - pm.KroneckerNormal, std_args, valuedomain=dom, - ref_rand=ref_rand_chol, extra_args=chol_args, - model_args=chol_args) + pm.KroneckerNormal, + std_args, + valuedomain=dom, + ref_rand=ref_rand_chol, + extra_args=chol_args, + model_args=chol_args, + ) pymc3_random( - pm.KroneckerNormal, std_args, valuedomain=dom, - ref_rand=ref_rand_evd, extra_args=evd_args, - model_args=evd_args) + pm.KroneckerNormal, + std_args, + valuedomain=dom, + ref_rand=ref_rand_evd, + extra_args=evd_args, + model_args=evd_args, + ) def test_mv_t(self): def ref_rand(size, nu, Sigma, mu): normal = st.multivariate_normal.rvs(cov=Sigma, size=size).T chi2 = st.chi2.rvs(df=nu, size=size) return mu + np.sqrt(nu) * (normal / chi2).T + for n in [2, 3]: - pymc3_random(pm.MvStudentT, - {'nu': Domain([5, 10, 25, 50]), 'Sigma': PdMatrix( - n), 'mu': Vector(R, n)}, - size=100, valuedomain=Vector(R, n), ref_rand=ref_rand) + pymc3_random( + pm.MvStudentT, + { + "nu": Domain([5, 10, 25, 50]), + "Sigma": PdMatrix(n), + "mu": Vector(R, n), + }, + size=100, + valuedomain=Vector(R, n), + ref_rand=ref_rand, + ) def test_dirichlet(self): def ref_rand(size, a): return st.dirichlet.rvs(a, size=size) + for n in [2, 3]: - pymc3_random(pm.Dirichlet, {'a': Vector(Rplus, n)}, - valuedomain=Simplex(n), size=100, ref_rand=ref_rand) + pymc3_random( + pm.Dirichlet, + {"a": Vector(Rplus, n)}, + valuedomain=Simplex(n), + size=100, + ref_rand=ref_rand, + ) def test_multinomial(self): def ref_rand(size, p, n): return nr.multinomial(pvals=p, n=n, size=size) + for n in [2, 3]: - pymc3_random_discrete(pm.Multinomial, {'p': Simplex(n), 'n': Nat}, - valuedomain=Vector(Nat, n), size=100, ref_rand=ref_rand) + pymc3_random_discrete( + pm.Multinomial, + {"p": Simplex(n), "n": Nat}, + valuedomain=Vector(Nat, n), + size=100, + ref_rand=ref_rand, + ) def test_gumbel(self): def ref_rand(size, mu, beta): return st.gumbel_r.rvs(loc=mu, scale=beta, size=size) - pymc3_random(pm.Gumbel, {'mu': R, 'beta': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.Gumbel, {"mu": R, "beta": Rplus}, ref_rand=ref_rand) def test_logistic(self): def ref_rand(size, mu, s): return st.logistic.rvs(loc=mu, scale=s, size=size) - pymc3_random(pm.Logistic, {'mu': R, 's': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.Logistic, {"mu": R, "s": Rplus}, ref_rand=ref_rand) def test_logitnormal(self): def ref_rand(size, mu, sigma): return expit(st.norm.rvs(loc=mu, scale=sigma, size=size)) - pymc3_random(pm.LogitNormal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand) + + pymc3_random(pm.LogitNormal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand) def test_moyal(self): def ref_rand(size, mu, sigma): return st.moyal.rvs(loc=mu, scale=sigma, size=size) - pymc3_random(pm.Moyal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand) - - @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32") + pymc3_random(pm.Moyal, {"mu": R, "sigma": Rplus}, ref_rand=ref_rand) + + @pytest.mark.xfail( + condition=(theano.config.floatX == "float32"), reason="Fails on float32" + ) def test_interpolated(self): for mu in R.vals: for sigma in Rplus.vals: - #pylint: disable=cell-var-from-loop + # pylint: disable=cell-var-from-loop def ref_rand(size): return st.norm.rvs(loc=mu, scale=sigma, size=size) - class TestedInterpolated (pm.Interpolated): - + class TestedInterpolated(pm.Interpolated): def __init__(self, **kwargs): x_points = np.linspace(mu - 5 * sigma, mu + 5 * sigma, 100) pdf_points = st.norm.pdf(x_points, loc=mu, scale=sigma) super().__init__( - x_points=x_points, - pdf_points=pdf_points, - **kwargs + x_points=x_points, pdf_points=pdf_points, **kwargs ) pymc3_random(TestedInterpolated, {}, ref_rand=ref_rand) - @pytest.mark.skip('Wishart random sampling not implemented.\n' - 'See https://github.com/pymc-devs/pymc3/issues/538') + @pytest.mark.skip( + "Wishart random sampling not implemented.\n" + "See https://github.com/pymc-devs/pymc3/issues/538" + ) def test_wishart(self): # Wishart non current recommended for use: # https://github.com/pymc-devs/pymc3/issues/538 @@ -878,85 +1092,77 @@ def test_wishart(self): def test_lkj(self): for n in [2, 10, 50]: - #pylint: disable=cell-var-from-loop - shape = n*(n-1)//2 + # pylint: disable=cell-var-from-loop + shape = n * (n - 1) // 2 def ref_rand(size, eta): - beta = eta - 1 + n/2 - return (st.beta.rvs(size=(size, shape), a=beta, b=beta)-.5)*2 - - class TestedLKJCorr (pm.LKJCorr): + beta = eta - 1 + n / 2 + return (st.beta.rvs(size=(size, shape), a=beta, b=beta) - 0.5) * 2 + class TestedLKJCorr(pm.LKJCorr): def __init__(self, **kwargs): - kwargs.pop('shape', None) + kwargs.pop("shape", None) super().__init__(n=n, **kwargs) - pymc3_random(TestedLKJCorr, - {'eta': Domain([1., 10., 100.])}, - size=10000//n, - ref_rand=ref_rand) + pymc3_random( + TestedLKJCorr, + {"eta": Domain([1.0, 10.0, 100.0])}, + size=10000 // n, + ref_rand=ref_rand, + ) def test_normalmixture(self): def ref_rand(size, w, mu, sigma): component = np.random.choice(w.size, size=size, p=w) return np.random.normal(mu[component], sigma[component], size=size) - pymc3_random(pm.NormalMixture, {'w': Simplex(2), - 'mu': Domain([[.05, 2.5], [-5., 1.]], edges=(None, None)), - 'sigma': Domain([[1, 1], [1.5, 2.]], edges=(None, None))}, - extra_args={'comp_shape': 2}, - size=1000, - ref_rand=ref_rand) - pymc3_random(pm.NormalMixture, {'w': Simplex(3), - 'mu': Domain([[-5., 1., 2.5]], edges=(None, None)), - 'sigma': Domain([[1.5, 2., 3.]], edges=(None, None))}, - extra_args={'comp_shape': 3}, - size=1000, - ref_rand=ref_rand) + pymc3_random( + pm.NormalMixture, + { + "w": Simplex(2), + "mu": Domain([[0.05, 2.5], [-5.0, 1.0]], edges=(None, None)), + "sigma": Domain([[1, 1], [1.5, 2.0]], edges=(None, None)), + }, + extra_args={"comp_shape": 2}, + size=1000, + ref_rand=ref_rand, + ) + pymc3_random( + pm.NormalMixture, + { + "w": Simplex(3), + "mu": Domain([[-5.0, 1.0, 2.5]], edges=(None, None)), + "sigma": Domain([[1.5, 2.0, 3.0]], edges=(None, None)), + }, + extra_args={"comp_shape": 3}, + size=1000, + ref_rand=ref_rand, + ) def test_mixture_random_shape(): # test the shape broadcasting in mixture random - y = np.concatenate([nr.poisson(5, size=10), - nr.poisson(9, size=10)]) + y = np.concatenate([nr.poisson(5, size=10), nr.poisson(9, size=10)]) with pm.Model() as m: comp0 = pm.Poisson.dist(mu=np.ones(2)) - w0 = pm.Dirichlet('w0', a=np.ones(2), shape=(2,)) - like0 = pm.Mixture('like0', - w=w0, - comp_dists=comp0, - observed=y) - - comp1 = pm.Poisson.dist(mu=np.ones((20, 2)), - shape=(20, 2)) - w1 = pm.Dirichlet('w1', a=np.ones(2), shape=(2,)) - like1 = pm.Mixture('like1', - w=w1, - comp_dists=comp1, - observed=y) + w0 = pm.Dirichlet("w0", a=np.ones(2), shape=(2,)) + like0 = pm.Mixture("like0", w=w0, comp_dists=comp0, observed=y) + + comp1 = pm.Poisson.dist(mu=np.ones((20, 2)), shape=(20, 2)) + w1 = pm.Dirichlet("w1", a=np.ones(2), shape=(2,)) + like1 = pm.Mixture("like1", w=w1, comp_dists=comp1, observed=y) comp2 = pm.Poisson.dist(mu=np.ones(2)) - w2 = pm.Dirichlet('w2', - a=np.ones(2), - shape=(20, 2)) - like2 = pm.Mixture('like2', - w=w2, - comp_dists=comp2, - observed=y) - - comp3 = pm.Poisson.dist(mu=np.ones(2), - shape=(20, 2)) - w3 = pm.Dirichlet('w3', - a=np.ones(2), - shape=(20, 2)) - like3 = pm.Mixture('like3', - w=w3, - comp_dists=comp3, - observed=y) - - rand0, rand1, rand2, rand3 = draw_values([like0, like1, like2, like3], - point=m.test_point, - size=100) + w2 = pm.Dirichlet("w2", a=np.ones(2), shape=(20, 2)) + like2 = pm.Mixture("like2", w=w2, comp_dists=comp2, observed=y) + + comp3 = pm.Poisson.dist(mu=np.ones(2), shape=(20, 2)) + w3 = pm.Dirichlet("w3", a=np.ones(2), shape=(20, 2)) + like3 = pm.Mixture("like3", w=w3, comp_dists=comp3, observed=y) + + rand0, rand1, rand2, rand3 = draw_values( + [like0, like1, like2, like3], point=m.test_point, size=100 + ) assert rand0.shape == (100, 20) assert rand1.shape == (100, 20) assert rand2.shape == (100, 20) @@ -964,54 +1170,36 @@ def test_mixture_random_shape(): with m: ppc = pm.sample_posterior_predictive([m.test_point], samples=200) - assert ppc['like0'].shape == (200, 20) - assert ppc['like1'].shape == (200, 20) - assert ppc['like2'].shape == (200, 20) - assert ppc['like3'].shape == (200, 20) + assert ppc["like0"].shape == (200, 20) + assert ppc["like1"].shape == (200, 20) + assert ppc["like2"].shape == (200, 20) + assert ppc["like3"].shape == (200, 20) + @pytest.mark.xfail def test_mixture_random_shape_fast(): # test the shape broadcasting in mixture random - y = np.concatenate([nr.poisson(5, size=10), - nr.poisson(9, size=10)]) + y = np.concatenate([nr.poisson(5, size=10), nr.poisson(9, size=10)]) with pm.Model() as m: comp0 = pm.Poisson.dist(mu=np.ones(2)) - w0 = pm.Dirichlet('w0', a=np.ones(2), shape=(2,)) - like0 = pm.Mixture('like0', - w=w0, - comp_dists=comp0, - observed=y) - - comp1 = pm.Poisson.dist(mu=np.ones((20, 2)), - shape=(20, 2)) - w1 = pm.Dirichlet('w1', a=np.ones(2), shape=(2,)) - like1 = pm.Mixture('like1', - w=w1, - comp_dists=comp1, - observed=y) + w0 = pm.Dirichlet("w0", a=np.ones(2), shape=(2,)) + like0 = pm.Mixture("like0", w=w0, comp_dists=comp0, observed=y) + + comp1 = pm.Poisson.dist(mu=np.ones((20, 2)), shape=(20, 2)) + w1 = pm.Dirichlet("w1", a=np.ones(2), shape=(2,)) + like1 = pm.Mixture("like1", w=w1, comp_dists=comp1, observed=y) comp2 = pm.Poisson.dist(mu=np.ones(2)) - w2 = pm.Dirichlet('w2', - a=np.ones(2), - shape=(20, 2)) - like2 = pm.Mixture('like2', - w=w2, - comp_dists=comp2, - observed=y) - - comp3 = pm.Poisson.dist(mu=np.ones(2), - shape=(20, 2)) - w3 = pm.Dirichlet('w3', - a=np.ones(2), - shape=(20, 2)) - like3 = pm.Mixture('like3', - w=w3, - comp_dists=comp3, - observed=y) - - rand0, rand1, rand2, rand3 = draw_values([like0, like1, like2, like3], - point=m.test_point, - size=100) + w2 = pm.Dirichlet("w2", a=np.ones(2), shape=(20, 2)) + like2 = pm.Mixture("like2", w=w2, comp_dists=comp2, observed=y) + + comp3 = pm.Poisson.dist(mu=np.ones(2), shape=(20, 2)) + w3 = pm.Dirichlet("w3", a=np.ones(2), shape=(20, 2)) + like3 = pm.Mixture("like3", w=w3, comp_dists=comp3, observed=y) + + rand0, rand1, rand2, rand3 = draw_values( + [like0, like1, like2, like3], point=m.test_point, size=100 + ) assert rand0.shape == (100, 20) assert rand1.shape == (100, 20) assert rand2.shape == (100, 20) @@ -1021,139 +1209,151 @@ def test_mixture_random_shape_fast(): # but I could be wrong. [2019/08/22:rpg] with m: ppc = pm.fast_sample_posterior_predictive([m.test_point], samples=200) - assert ppc['like0'].shape == (200, 20) - assert ppc['like1'].shape == (200, 20) - assert ppc['like2'].shape == (200, 20) - assert ppc['like3'].shape == (200, 20) - + assert ppc["like0"].shape == (200, 20) + assert ppc["like1"].shape == (200, 20) + assert ppc["like2"].shape == (200, 20) + assert ppc["like3"].shape == (200, 20) -class TestDensityDist(): +class TestDensityDist: @pytest.mark.parametrize("shape", [(), (3,), (3, 2)], ids=str) def test_density_dist_with_random_sampleable(self, shape): with pm.Model() as model: - mu = pm.Normal('mu', 0, 1) + mu = pm.Normal("mu", 0, 1) normal_dist = pm.Normal.dist(mu, 1, shape=shape) obs = pm.DensityDist( - 'density_dist', + "density_dist", normal_dist.logp, observed=np.random.randn(100, *shape), shape=shape, - random=normal_dist.random) + random=normal_dist.random, + ) trace = pm.sample(100) samples = 500 size = 100 - ppc = pm.sample_posterior_predictive(trace, samples=samples, model=model, size=size) - assert ppc['density_dist'].shape == (samples, size) + obs.distribution.shape + ppc = pm.sample_posterior_predictive( + trace, samples=samples, model=model, size=size + ) + assert ppc["density_dist"].shape == (samples, size) + obs.distribution.shape # ppc = pm.fast_sample_posterior_predictive(trace, samples=samples, model=model, size=size) # assert ppc['density_dist'].shape == (samples, size) + obs.distribution.shape - @pytest.mark.parametrize("shape", [(), (3,), (3, 2)], ids=str) def test_density_dist_with_random_sampleable_failure(self, shape): with pm.Model() as model: - mu = pm.Normal('mu', 0, 1) + mu = pm.Normal("mu", 0, 1) normal_dist = pm.Normal.dist(mu, 1, shape=shape) pm.DensityDist( - 'density_dist', + "density_dist", normal_dist.logp, observed=np.random.randn(100, *shape), shape=shape, random=normal_dist.random, - wrap_random_with_dist_shape=False + wrap_random_with_dist_shape=False, ) trace = pm.sample(100) samples = 500 with pytest.raises(RuntimeError): - pm.sample_posterior_predictive(trace, samples=samples, model=model, size=100) + pm.sample_posterior_predictive( + trace, samples=samples, model=model, size=100 + ) with pytest.raises((TypeError, RuntimeError)): - pm.fast_sample_posterior_predictive(trace, samples=samples, model=model, size=100) - + pm.fast_sample_posterior_predictive( + trace, samples=samples, model=model, size=100 + ) @pytest.mark.parametrize("shape", [(), (3,), (3, 2)], ids=str) def test_density_dist_with_random_sampleable_hidden_error(self, shape): with pm.Model() as model: - mu = pm.Normal('mu', 0, 1) + mu = pm.Normal("mu", 0, 1) normal_dist = pm.Normal.dist(mu, 1, shape=shape) obs = pm.DensityDist( - 'density_dist', + "density_dist", normal_dist.logp, observed=np.random.randn(100, *shape), shape=shape, random=normal_dist.random, wrap_random_with_dist_shape=False, - check_shape_in_random=False + check_shape_in_random=False, ) trace = pm.sample(100) samples = 500 ppc = pm.sample_posterior_predictive(trace, samples=samples, model=model) - assert len(ppc['density_dist']) == samples - assert ((samples,) + obs.distribution.shape) != ppc['density_dist'].shape + assert len(ppc["density_dist"]) == samples + assert ((samples,) + obs.distribution.shape) != ppc["density_dist"].shape ppc = pm.fast_sample_posterior_predictive(trace, samples=samples, model=model) - assert len(ppc['density_dist']) == samples - assert ((samples,) + obs.distribution.shape) != ppc['density_dist'].shape - + assert len(ppc["density_dist"]) == samples + assert ((samples,) + obs.distribution.shape) != ppc["density_dist"].shape def test_density_dist_with_random_sampleable_handcrafted_success(self): with pm.Model() as model: - mu = pm.Normal('mu', 0, 1) + mu = pm.Normal("mu", 0, 1) normal_dist = pm.Normal.dist(mu, 1) rvs = pm.Normal.dist(mu, 1, shape=100).random obs = pm.DensityDist( - 'density_dist', + "density_dist", normal_dist.logp, observed=np.random.randn(100), random=rvs, - wrap_random_with_dist_shape=False + wrap_random_with_dist_shape=False, ) trace = pm.sample(100) samples = 500 size = 100 - ppc = pm.sample_posterior_predictive(trace, samples=samples, model=model, size=size) - assert ppc['density_dist'].shape == (samples, size) + obs.distribution.shape + ppc = pm.sample_posterior_predictive( + trace, samples=samples, model=model, size=size + ) + assert ppc["density_dist"].shape == (samples, size) + obs.distribution.shape @pytest.mark.xfail def test_density_dist_with_random_sampleable_handcrafted_success_fast(self): with pm.Model() as model: - mu = pm.Normal('mu', 0, 1) + mu = pm.Normal("mu", 0, 1) normal_dist = pm.Normal.dist(mu, 1) rvs = pm.Normal.dist(mu, 1, shape=100).random obs = pm.DensityDist( - 'density_dist', + "density_dist", normal_dist.logp, observed=np.random.randn(100), random=rvs, - wrap_random_with_dist_shape=False + wrap_random_with_dist_shape=False, ) trace = pm.sample(100) samples = 500 size = 100 - ppc = pm.fast_sample_posterior_predictive(trace, samples=samples, model=model, size=size) - assert ppc['density_dist'].shape == (samples, size) + obs.distribution.shape - + ppc = pm.fast_sample_posterior_predictive( + trace, samples=samples, model=model, size=size + ) + assert ppc["density_dist"].shape == (samples, size) + obs.distribution.shape def test_density_dist_without_random_not_sampleable(self): with pm.Model() as model: - mu = pm.Normal('mu', 0, 1) + mu = pm.Normal("mu", 0, 1) normal_dist = pm.Normal.dist(mu, 1) - pm.DensityDist('density_dist', normal_dist.logp, observed=np.random.randn(100)) + pm.DensityDist( + "density_dist", normal_dist.logp, observed=np.random.randn(100) + ) trace = pm.sample(100) samples = 500 with pytest.raises(ValueError): - pm.sample_posterior_predictive(trace, samples=samples, model=model, size=100) + pm.sample_posterior_predictive( + trace, samples=samples, model=model, size=100 + ) with pytest.raises((TypeError, ValueError)): - pm.fast_sample_posterior_predictive(trace, samples=samples, model=model, size=100) + pm.fast_sample_posterior_predictive( + trace, samples=samples, model=model, size=100 + ) class TestNestedRandom(SeededTest): @@ -1163,7 +1363,7 @@ def build_model(self, distribution, shape, nested_rvs_info): for rv_name, info in nested_rvs_info.items(): try: value, nested_shape = info - loc = 0. + loc = 0.0 except ValueError: value, nested_shape, loc = info if value is None: @@ -1182,13 +1382,7 @@ def build_model(self, distribution, shape, nested_rvs_info): ) return model, rv, nested_rvs - def sample_prior( - self, - distribution, - shape, - nested_rvs_info, - prior_samples - ): + def sample_prior(self, distribution, shape, nested_rvs_info, prior_samples): model, rv, nested_rvs = self.build_model( distribution, shape, @@ -1202,8 +1396,24 @@ def sample_prior( [ [10, (3,), (None, tuple()), (None, (3,))], [10, (3,), (None, (3,)), (None, tuple())], - [10, (4, 3,), (None, (3,)), (None, (3,))], - [10, (4, 3,), (None, (3,)), (None, (4, 3))], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (None, (3,)), + ], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (None, (4, 3)), + ], ], ids=str, ) @@ -1229,8 +1439,26 @@ def test_NegativeBinomial( [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))], + [ + 10, + ( + 4, + 3, + ), + (0.5, (3,)), + (None, (3,)), + (None, (3,)), + ], + [ + 10, + ( + 4, + 3, + ), + (0.5, (3,)), + (None, (3,)), + (None, (4, 3)), + ], ], ids=str, ) @@ -1257,8 +1485,24 @@ def test_ZeroInflatedNegativeBinomial( [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))], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (None, (3,)), + ], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (None, (4, 3)), + ], ], ids=str, ) @@ -1280,18 +1524,114 @@ def test_Rice( @pytest.mark.parametrize( ["prior_samples", "shape", "mu", "sigma", "lower", "upper"], [ - [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))], + [ + 10, + (3,), + (None, tuple()), + (1.0, tuple()), + (None, tuple(), -1), + (None, (3,)), + ], + [ + 10, + (3,), + (None, tuple()), + (1.0, tuple()), + (None, tuple(), -1), + (None, (3,)), + ], + [ + 10, + (3,), + (None, tuple()), + (1.0, tuple()), + (None, (3,), -1), + (None, tuple()), + ], + [ + 10, + (3,), + (None, tuple()), + (1.0, tuple()), + (None, (3,), -1), + (None, tuple()), + ], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (1.0, tuple()), + (None, (3,), -1), + (None, (3,)), + ], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (1.0, tuple()), + (None, (3,), -1), + (None, (4, 3)), + ], + [ + 10, + (3,), + (0.0, tuple()), + (None, tuple()), + (None, tuple(), -1), + (None, (3,)), + ], + [ + 10, + (3,), + (0.0, tuple()), + (None, tuple()), + (None, tuple(), -1), + (None, (3,)), + ], + [ + 10, + (3,), + (0.0, tuple()), + (None, tuple()), + (None, (3,), -1), + (None, tuple()), + ], + [ + 10, + (3,), + (0.0, tuple()), + (None, tuple()), + (None, (3,), -1), + (None, tuple()), + ], + [ + 10, + ( + 4, + 3, + ), + (0.0, tuple()), + (None, (3,)), + (None, (3,), -1), + (None, (3,)), + ], + [ + 10, + ( + 4, + 3, + ), + (0.0, tuple()), + (None, (3,)), + (None, (3,), -1), + (None, (4, 3)), + ], ], ids=str, ) @@ -1312,15 +1652,32 @@ def test_TruncatedNormal( ) 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)], + [10, (3,), (None, tuple()), (-1.0, (3,)), (2, tuple())], + [10, (3,), (None, tuple()), (-1.0, tuple()), (None, tuple(), 1)], + [10, (3,), (None, (3,)), (-1.0, tuple()), (None, tuple(), 1)], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (-1.0, tuple()), + (None, (3,), 1), + ], + [ + 10, + ( + 4, + 3, + ), + (None, (3,)), + (None, tuple(), -1), + (None, (3,), 1), + ], ], ids=str, ) From d797ef7d3a6c55a264896a26122436bff5392839 Mon Sep 17 00:00:00 2001 From: harivallabha Date: Fri, 18 Sep 2020 01:43:21 +0530 Subject: [PATCH 7/8] Remove unnecessary bound condition in logprob of hypergeom --- pymc3/distributions/discrete.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index 919d75c059..68e545d407 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -948,7 +948,6 @@ def logp(self, value): 0 <= k, k <= N, 0 <= n, - 0 <= N, n - N + k <= value, 0 <= value, value <= k, From c708e2fcabed3eb676beb5f22aae19b6ccbb11d0 Mon Sep 17 00:00:00 2001 From: harivallabha Date: Fri, 18 Sep 2020 01:59:08 +0530 Subject: [PATCH 8/8] Change logp of hypergeometric to mirror the scipy implementation --- pymc3/distributions/discrete.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index 68e545d407..da2023999d 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -943,16 +943,17 @@ def logp(self, value): N = self.N k = self.k n = self.n - return bound( - binomln(k, value) + binomln(N - k, n - value) - binomln(N, n), - 0 <= k, - k <= N, - 0 <= n, - n - N + k <= value, - 0 <= value, - value <= k, - value <= n, + tot, good = N, k + bad = tot - good + result = ( + betaln(good + 1, 1) + + betaln(bad + 1, 1) + + betaln(tot - n + 1, n + 1) + - betaln(value + 1, good - value + 1) + - betaln(n - value + 1, bad - n + value + 1) + - betaln(tot + 1, 1) ) + return result def _repr_latex_(self, name=None, dist=None): if dist is None: