From f1a67b90373592e331cdce23d93e1276dce0db73 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 2 Feb 2022 16:47:02 +1100 Subject: [PATCH 01/15] Initial commit --- pymc_experimental/distributions/__init__.py | 26 ++- pymc_experimental/distributions/genextreme.py | 217 ++++++++++++++++++ pymc_experimental/tests/test_distributions.py | 61 +++++ .../tests/test_distributions_random.py | 40 ++++ 4 files changed, 342 insertions(+), 2 deletions(-) create mode 100644 pymc_experimental/distributions/genextreme.py create mode 100644 pymc_experimental/tests/test_distributions.py create mode 100644 pymc_experimental/tests/test_distributions_random.py diff --git a/pymc_experimental/distributions/__init__.py b/pymc_experimental/distributions/__init__.py index 41ff5483..3eb86c72 100644 --- a/pymc_experimental/distributions/__init__.py +++ b/pymc_experimental/distributions/__init__.py @@ -1,2 +1,24 @@ -from pymc_experimental.distributions import histogram_utils -from pymc_experimental.distributions.histogram_utils import histogram_approximation +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# coding: utf-8 +""" +Experimental probability distributions for stochastic nodes in PyMC. +""" + +from pymc-experimental.distributions.genextreme import GenExtreme + +__all__ = [ + "GenExtreme", +] diff --git a/pymc_experimental/distributions/genextreme.py b/pymc_experimental/distributions/genextreme.py new file mode 100644 index 00000000..da0b7b93 --- /dev/null +++ b/pymc_experimental/distributions/genextreme.py @@ -0,0 +1,217 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# coding: utf-8 +""" +Experimental probability distributions for stochastic nodes in PyMC. +""" + +from typing import List, Tuple + +import aesara +import aesara.tensor as at +import numpy as np +from scipy import stats + +from aesara.tensor.random.op import RandomVariable +from pymc.distributions.distribution import Continuous +from aesara.tensor.var import TensorVariable +from pymc.aesaraf import floatX +from pymc.distributions.dist_math import bound + + +class GenExtremeRV(RandomVariable): + name: str = "Generalized Extreme Value" + ndim_supp: int = 0 + ndims_params: List[int] = [0, 0, 0] + dtype: str = "floatX" + _print_name: Tuple[str, str] = ("Generalized Extreme Value", "\\operatorname{GEV}") + + def __call__( + self, mu=0.0, sigma=1.0, xi=0.0, size=None, **kwargs + ) -> TensorVariable: + return super().__call__(mu, sigma, xi, size=size, **kwargs) + + @classmethod + def rng_fn( + cls, + rng: np.random.RandomState, + mu: np.ndarray, + sigma: np.ndarray, + xi: np.ndarray, + size: Tuple[int, ...], + ) -> np.ndarray: + # Notice negative here, since remainder of GenExtreme is based on Coles parametrization + return stats.genextreme.rvs( + c=-xi, loc=mu, scale=sigma, random_state=rng, size=size + ) + + +gev = GenExtremeRV() + + +class GenExtreme(Continuous): + r""" + Univariate Generalized Extreme Value log-likelihood + + The cdf of this distribution is + + .. math:: + + G(x \mid \mu, \sigma, \xi) = \exp\left[ -\left(1 + \xi z\right)^{-\frac{1}{\xi}} \right] + + where + + .. math:: + + z = \frac{x - \mu}{\sigma} + + and is defined on the set: + + .. math:: + + \left\{x: 1 + \xi\left(\frac{x-\mu}{\sigma}\right) > 0 \right\}. + + Note that this parametrization is per Coles (2001), and differs from that of + Scipy in the sign of the shape parameter, :math:`\xi`. + + .. plot:: + + import matplotlib.pyplot as plt + import numpy as np + import scipy.stats as st + import arviz as az + plt.style.use('arviz-darkgrid') + x = np.linspace(-10, 20, 200) + mus = [0., 4., -1.] + sigmas = [2., 2., 4.] + xis = [-0.3, 0.0, 0.3] + for mu, sigma, xi in zip(mus, sigmas, xis): + pdf = st.genextreme.pdf(x, c=-xi, loc=mu, scale=sigma) + plt.plot(x, pdf, label=rf'$\mu$ = {mu}, $\sigma$ = {sigma}, $\xi$={xi}') + plt.xlabel('x', fontsize=12) + plt.ylabel('f(x)', fontsize=12) + plt.legend(loc=1) + plt.show() + + + ======== ========================================================================= + Support * :math:`x \in [\mu - \sigma/\xi, +\infty]`, when :math:`\xi > 0` + * :math:`x \in \mathbb{R}` when :math:`\xi = 0` + * :math:`x \in [-\infty, \mu - \sigma/\xi]`, when :math:`\xi < 0` + Mean * :math:`\mu + \sigma(g_1 - 1)/\xi`, when :math:`\xi \neq 0, \xi < 1` + * :math:`\mu + \sigma \gamma`, when :math:`\xi = 0` + * :math:`\infty`, when :math:`\xi \geq 1` + where :math:`\gamma` is the Euler-Mascheroni constant, and + :math:`g_k = \Gamma (1-k\xi)` + Variance * :math:`\sigma^2 (g_2 - g_1^2)/\xi^2`, when :math:`\xi \neq 0, \xi < 0.5` + * :math:`\frac{\pi^2}{6} \sigma^2`, when :math:`\xi = 0` + * :math:`\infty`, when :math:`\xi \geq 0.5` + ======== ========================================================================= + + Parameters + ---------- + mu: float + Location parameter. + sigma: float + Scale parameter (sigma > 0). + xi: float + Shape parameter + scipy: bool + Whether or not to use the Scipy interpretation of the shape parameter + (defaults to `False`). + + References + ---------- + .. [Coles2001] Coles, S.G. (2001). + An Introduction to the Statistical Modeling of Extreme Values + Springer-Verlag, London + + """ + + rv_op = gev + + @classmethod + def dist(cls, mu=0, sigma=1, xi=0, scipy=False, **kwargs): + # If SciPy, use its parametrization, otherwise convert to standard + if scipy: + xi = -xi + mu = at.as_tensor_variable(floatX(mu)) + sigma = at.as_tensor_variable(floatX(sigma)) + xi = at.as_tensor_variable(floatX(xi)) + + return super().dist([mu, sigma, xi], **kwargs) + + def logp(value, mu, sigma, xi): + """ + Calculate log-probability of Generalized Extreme Value 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 Aesara tensor + + Returns + ------- + TensorVariable + """ + scaled = (value - mu) / sigma + + logp_expression = at.switch( + at.isclose(xi, 0), + -at.log(sigma) - scaled - at.exp(-scaled), + -at.log(sigma) + - ((xi + 1) / xi) * at.log1p(xi * scaled) + - at.pow(1 + xi * scaled, -1 / xi), + ) + # bnd = mu - sigma/xi + return bound( + logp_expression, + 1 + xi * (value - mu) / sigma > 0, + # at.switch(xi > 0, value > bnd, value < bnd), + sigma > 0, + ) + + def logcdf(value, mu, sigma, xi): + """ + Compute the log of the cumulative distribution function for Generalized Extreme Value + distribution at the specified value. + + Parameters + ---------- + value: numeric or np.ndarray or `TensorVariable` + Value(s) for which log CDF is calculated. If the log CDF for + multiple values are desired the values must be provided in a numpy + array or `TensorVariable`. + + Returns + ------- + TensorVariable + """ + scaled = (value - mu) / sigma + logc_expression = at.switch( + at.isclose(xi, 0), -at.exp(-scaled), -at.pow(1 + xi * scaled, -1 / xi) + ) + return bound(logc_expression, 1 + xi * scaled > 0, sigma > 0) + + def get_moment(value_var, size, mu, sigma, xi): + r""" + Using the mode, as the mean can be infinite when :math:`\xi > 1` + """ + mode = at.switch( + at.isclose(xi, 0), mu, mu + sigma * (at.pow(1 + xi, -xi) - 1) / xi + ) + return at.full(size, mode, dtype=aesara.config.floatX) diff --git a/pymc_experimental/tests/test_distributions.py b/pymc_experimental/tests/test_distributions.py new file mode 100644 index 00000000..479f9a18 --- /dev/null +++ b/pymc_experimental/tests/test_distributions.py @@ -0,0 +1,61 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import functools +import itertools +import sys + +import aesara +import aesara.tensor as at +import numpy as np +import numpy.random as nr + +import pytest +import scipy.stats +import scipy.stats.distributions as sp + +from aesara.compile.mode import Mode +from aesara.graph.basic import ancestors +from aesara.tensor.random.op import RandomVariable +from aesara.tensor.var import TensorVariable +from numpy import array, inf, log +from numpy.testing import assert_allclose, assert_almost_equal, assert_equal +from scipy import integrate +from scipy.special import erf, logit + +import pymc as pm + +from pymc.aesaraf import floatX, intX +from pymc-experimental.distributions import ( + GenExtreme, +) +from pymc.distributions.shape_utils import to_tuple +from pymc.math import kronecker +from pymc.model import Deterministic, Model, Point, Potential +from pymc.tests.helpers import select_by_precision +from pymc.vartypes import continuous_types + +def test_genextreme(self): + self.check_logp( + GenExtreme, + R, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), + ) + self.check_logcdf( + GenExtreme, + R, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), + ) + diff --git a/pymc_experimental/tests/test_distributions_random.py b/pymc_experimental/tests/test_distributions_random.py new file mode 100644 index 00000000..76a5dfbb --- /dev/null +++ b/pymc_experimental/tests/test_distributions_random.py @@ -0,0 +1,40 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import functools +import itertools + +from typing import Callable, List, Optional + +import aesara +import numpy as np +import numpy.random as nr +import numpy.testing as npt +import pytest +import scipy.stats as st + +from numpy.testing import assert_almost_equal, assert_array_almost_equal + + +class TestGenExtreme(BaseTestDistribution): + pymc_dist = pm.GenExtreme + pymc_dist_params = {"mu": 0, "sigma": 1, "xi": -0.1} + expected_rv_op_params = {"mu": 0, "sigma": 1, "xi": -0.1} + # Notice, using different parametrization of xi sign to scipy + reference_dist_params = {"loc": 0, "scale": 1, "c": 0.1} + reference_dist = seeded_scipy_distribution_builder("genextreme") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] From 959677f2f26ae778359556917aa033bb32c17041 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Fri, 4 Feb 2022 22:00:26 +1100 Subject: [PATCH 02/15] Updated structure; code operational --- pymc_experimental/distributions/continuous.py | 218 ++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 pymc_experimental/distributions/continuous.py diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py new file mode 100644 index 00000000..0b8e5246 --- /dev/null +++ b/pymc_experimental/distributions/continuous.py @@ -0,0 +1,218 @@ +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# coding: utf-8 +""" +Experimental probability distributions for stochastic nodes in PyMC. + +The imports from pymc are not fully replicated here: add imports as necessary. +""" + +from typing import List, Tuple +import aesara +import aesara.tensor as at +import numpy as np +from scipy import stats + +from aesara.tensor.random.op import RandomVariable +from pymc.distributions.distribution import Continuous +from aesara.tensor.var import TensorVariable +from pymc.aesaraf import floatX +from pymc.distributions.dist_math import check_parameters + + +class GenExtremeRV(RandomVariable): + name: str = "Generalized Extreme Value" + ndim_supp: int = 0 + ndims_params: List[int] = [0, 0, 0] + dtype: str = "floatX" + _print_name: Tuple[str, str] = ("Generalized Extreme Value", "\\operatorname{GEV}") + + def __call__( + self, mu=0.0, sigma=1.0, xi=0.0, size=None, **kwargs + ) -> TensorVariable: + return super().__call__(mu, sigma, xi, size=size, **kwargs) + + @classmethod + def rng_fn( + cls, + rng: np.random.RandomState, + mu: np.ndarray, + sigma: np.ndarray, + xi: np.ndarray, + size: Tuple[int, ...], + ) -> np.ndarray: + # Notice negative here, since remainder of GenExtreme is based on Coles parametrization + return stats.genextreme.rvs( + c=-xi, loc=mu, scale=sigma, random_state=rng, size=size + ) + + +gev = GenExtremeRV() + + +class GenExtreme(Continuous): + r""" + Univariate Generalized Extreme Value log-likelihood + + The cdf of this distribution is + + .. math:: + + G(x \mid \mu, \sigma, \xi) = \exp\left[ -\left(1 + \xi z\right)^{-\frac{1}{\xi}} \right] + + where + + .. math:: + + z = \frac{x - \mu}{\sigma} + + and is defined on the set: + + .. math:: + + \left\{x: 1 + \xi\left(\frac{x-\mu}{\sigma}\right) > 0 \right\}. + + Note that this parametrization is per Coles (2001), and differs from that of + Scipy in the sign of the shape parameter, :math:`\xi`. + + .. plot:: + + import matplotlib.pyplot as plt + import numpy as np + import scipy.stats as st + import arviz as az + plt.style.use('arviz-darkgrid') + x = np.linspace(-10, 20, 200) + mus = [0., 4., -1.] + sigmas = [2., 2., 4.] + xis = [-0.3, 0.0, 0.3] + for mu, sigma, xi in zip(mus, sigmas, xis): + pdf = st.genextreme.pdf(x, c=-xi, loc=mu, scale=sigma) + plt.plot(x, pdf, label=rf'$\mu$ = {mu}, $\sigma$ = {sigma}, $\xi$={xi}') + plt.xlabel('x', fontsize=12) + plt.ylabel('f(x)', fontsize=12) + plt.legend(loc=1) + plt.show() + + + ======== ========================================================================= + Support * :math:`x \in [\mu - \sigma/\xi, +\infty]`, when :math:`\xi > 0` + * :math:`x \in \mathbb{R}` when :math:`\xi = 0` + * :math:`x \in [-\infty, \mu - \sigma/\xi]`, when :math:`\xi < 0` + Mean * :math:`\mu + \sigma(g_1 - 1)/\xi`, when :math:`\xi \neq 0, \xi < 1` + * :math:`\mu + \sigma \gamma`, when :math:`\xi = 0` + * :math:`\infty`, when :math:`\xi \geq 1` + where :math:`\gamma` is the Euler-Mascheroni constant, and + :math:`g_k = \Gamma (1-k\xi)` + Variance * :math:`\sigma^2 (g_2 - g_1^2)/\xi^2`, when :math:`\xi \neq 0, \xi < 0.5` + * :math:`\frac{\pi^2}{6} \sigma^2`, when :math:`\xi = 0` + * :math:`\infty`, when :math:`\xi \geq 0.5` + ======== ========================================================================= + + Parameters + ---------- + mu: float + Location parameter. + sigma: float + Scale parameter (sigma > 0). + xi: float + Shape parameter + scipy: bool + Whether or not to use the Scipy interpretation of the shape parameter + (defaults to `False`). + + References + ---------- + .. [Coles2001] Coles, S.G. (2001). + An Introduction to the Statistical Modeling of Extreme Values + Springer-Verlag, London + + """ + + rv_op = gev + + @classmethod + def dist(cls, mu=0, sigma=1, xi=0, scipy=False, **kwargs): + # If SciPy, use its parametrization, otherwise convert to standard + if scipy: + xi = -xi + mu = at.as_tensor_variable(floatX(mu)) + sigma = at.as_tensor_variable(floatX(sigma)) + xi = at.as_tensor_variable(floatX(xi)) + + return super().dist([mu, sigma, xi], **kwargs) + + def logp(value, mu, sigma, xi): + """ + Calculate log-probability of Generalized Extreme Value 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 Aesara tensor + + Returns + ------- + TensorVariable + """ + scaled = (value - mu) / sigma + + logp_expression = at.switch( + at.isclose(xi, 0), + -at.log(sigma) - scaled - at.exp(-scaled), + -at.log(sigma) + - ((xi + 1) / xi) * at.log1p(xi * scaled) + - at.pow(1 + xi * scaled, -1 / xi), + ) + # bnd = mu - sigma/xi + return check_parameters( + logp_expression, + 1 + xi * (value - mu) / sigma > 0, + # at.switch(xi > 0, value > bnd, value < bnd), + sigma > 0, + ) + + def logcdf(value, mu, sigma, xi): + """ + Compute the log of the cumulative distribution function for Generalized Extreme Value + distribution at the specified value. + + Parameters + ---------- + value: numeric or np.ndarray or `TensorVariable` + Value(s) for which log CDF is calculated. If the log CDF for + multiple values are desired the values must be provided in a numpy + array or `TensorVariable`. + + Returns + ------- + TensorVariable + """ + scaled = (value - mu) / sigma + logc_expression = at.switch( + at.isclose(xi, 0), -at.exp(-scaled), -at.pow(1 + xi * scaled, -1 / xi) + ) + return check_parameters(logc_expression, 1 + xi * scaled > 0, sigma > 0) + + def get_moment(value_var, size, mu, sigma, xi): + r""" + Using the mode, as the mean can be infinite when :math:`\xi > 1` + """ + mode = at.switch( + at.isclose(xi, 0), mu, mu + sigma * (at.pow(1 + xi, -xi) - 1) / xi + ) + return at.full(size, mode, dtype=aesara.config.floatX) From 929bdfac79438a742b3946bfc2711aef9875c249 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 16 Feb 2022 14:50:03 +1100 Subject: [PATCH 03/15] Change scipy.ststats.distribution import --- pymc_experimental/distributions/__init__.py | 4 +- pymc_experimental/distributions/genextreme.py | 217 ------------------ pymc_experimental/tests/test_distributions.py | 77 +++---- .../tests/test_distributions_random.py | 20 +- 4 files changed, 47 insertions(+), 271 deletions(-) delete mode 100644 pymc_experimental/distributions/genextreme.py diff --git a/pymc_experimental/distributions/__init__.py b/pymc_experimental/distributions/__init__.py index 3eb86c72..c04f53de 100644 --- a/pymc_experimental/distributions/__init__.py +++ b/pymc_experimental/distributions/__init__.py @@ -17,7 +17,9 @@ Experimental probability distributions for stochastic nodes in PyMC. """ -from pymc-experimental.distributions.genextreme import GenExtreme +from pymc_experimental.distributions.continuous import ( + GenExtreme, +) __all__ = [ "GenExtreme", diff --git a/pymc_experimental/distributions/genextreme.py b/pymc_experimental/distributions/genextreme.py deleted file mode 100644 index da0b7b93..00000000 --- a/pymc_experimental/distributions/genextreme.py +++ /dev/null @@ -1,217 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -# coding: utf-8 -""" -Experimental probability distributions for stochastic nodes in PyMC. -""" - -from typing import List, Tuple - -import aesara -import aesara.tensor as at -import numpy as np -from scipy import stats - -from aesara.tensor.random.op import RandomVariable -from pymc.distributions.distribution import Continuous -from aesara.tensor.var import TensorVariable -from pymc.aesaraf import floatX -from pymc.distributions.dist_math import bound - - -class GenExtremeRV(RandomVariable): - name: str = "Generalized Extreme Value" - ndim_supp: int = 0 - ndims_params: List[int] = [0, 0, 0] - dtype: str = "floatX" - _print_name: Tuple[str, str] = ("Generalized Extreme Value", "\\operatorname{GEV}") - - def __call__( - self, mu=0.0, sigma=1.0, xi=0.0, size=None, **kwargs - ) -> TensorVariable: - return super().__call__(mu, sigma, xi, size=size, **kwargs) - - @classmethod - def rng_fn( - cls, - rng: np.random.RandomState, - mu: np.ndarray, - sigma: np.ndarray, - xi: np.ndarray, - size: Tuple[int, ...], - ) -> np.ndarray: - # Notice negative here, since remainder of GenExtreme is based on Coles parametrization - return stats.genextreme.rvs( - c=-xi, loc=mu, scale=sigma, random_state=rng, size=size - ) - - -gev = GenExtremeRV() - - -class GenExtreme(Continuous): - r""" - Univariate Generalized Extreme Value log-likelihood - - The cdf of this distribution is - - .. math:: - - G(x \mid \mu, \sigma, \xi) = \exp\left[ -\left(1 + \xi z\right)^{-\frac{1}{\xi}} \right] - - where - - .. math:: - - z = \frac{x - \mu}{\sigma} - - and is defined on the set: - - .. math:: - - \left\{x: 1 + \xi\left(\frac{x-\mu}{\sigma}\right) > 0 \right\}. - - Note that this parametrization is per Coles (2001), and differs from that of - Scipy in the sign of the shape parameter, :math:`\xi`. - - .. plot:: - - import matplotlib.pyplot as plt - import numpy as np - import scipy.stats as st - import arviz as az - plt.style.use('arviz-darkgrid') - x = np.linspace(-10, 20, 200) - mus = [0., 4., -1.] - sigmas = [2., 2., 4.] - xis = [-0.3, 0.0, 0.3] - for mu, sigma, xi in zip(mus, sigmas, xis): - pdf = st.genextreme.pdf(x, c=-xi, loc=mu, scale=sigma) - plt.plot(x, pdf, label=rf'$\mu$ = {mu}, $\sigma$ = {sigma}, $\xi$={xi}') - plt.xlabel('x', fontsize=12) - plt.ylabel('f(x)', fontsize=12) - plt.legend(loc=1) - plt.show() - - - ======== ========================================================================= - Support * :math:`x \in [\mu - \sigma/\xi, +\infty]`, when :math:`\xi > 0` - * :math:`x \in \mathbb{R}` when :math:`\xi = 0` - * :math:`x \in [-\infty, \mu - \sigma/\xi]`, when :math:`\xi < 0` - Mean * :math:`\mu + \sigma(g_1 - 1)/\xi`, when :math:`\xi \neq 0, \xi < 1` - * :math:`\mu + \sigma \gamma`, when :math:`\xi = 0` - * :math:`\infty`, when :math:`\xi \geq 1` - where :math:`\gamma` is the Euler-Mascheroni constant, and - :math:`g_k = \Gamma (1-k\xi)` - Variance * :math:`\sigma^2 (g_2 - g_1^2)/\xi^2`, when :math:`\xi \neq 0, \xi < 0.5` - * :math:`\frac{\pi^2}{6} \sigma^2`, when :math:`\xi = 0` - * :math:`\infty`, when :math:`\xi \geq 0.5` - ======== ========================================================================= - - Parameters - ---------- - mu: float - Location parameter. - sigma: float - Scale parameter (sigma > 0). - xi: float - Shape parameter - scipy: bool - Whether or not to use the Scipy interpretation of the shape parameter - (defaults to `False`). - - References - ---------- - .. [Coles2001] Coles, S.G. (2001). - An Introduction to the Statistical Modeling of Extreme Values - Springer-Verlag, London - - """ - - rv_op = gev - - @classmethod - def dist(cls, mu=0, sigma=1, xi=0, scipy=False, **kwargs): - # If SciPy, use its parametrization, otherwise convert to standard - if scipy: - xi = -xi - mu = at.as_tensor_variable(floatX(mu)) - sigma = at.as_tensor_variable(floatX(sigma)) - xi = at.as_tensor_variable(floatX(xi)) - - return super().dist([mu, sigma, xi], **kwargs) - - def logp(value, mu, sigma, xi): - """ - Calculate log-probability of Generalized Extreme Value 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 Aesara tensor - - Returns - ------- - TensorVariable - """ - scaled = (value - mu) / sigma - - logp_expression = at.switch( - at.isclose(xi, 0), - -at.log(sigma) - scaled - at.exp(-scaled), - -at.log(sigma) - - ((xi + 1) / xi) * at.log1p(xi * scaled) - - at.pow(1 + xi * scaled, -1 / xi), - ) - # bnd = mu - sigma/xi - return bound( - logp_expression, - 1 + xi * (value - mu) / sigma > 0, - # at.switch(xi > 0, value > bnd, value < bnd), - sigma > 0, - ) - - def logcdf(value, mu, sigma, xi): - """ - Compute the log of the cumulative distribution function for Generalized Extreme Value - distribution at the specified value. - - Parameters - ---------- - value: numeric or np.ndarray or `TensorVariable` - Value(s) for which log CDF is calculated. If the log CDF for - multiple values are desired the values must be provided in a numpy - array or `TensorVariable`. - - Returns - ------- - TensorVariable - """ - scaled = (value - mu) / sigma - logc_expression = at.switch( - at.isclose(xi, 0), -at.exp(-scaled), -at.pow(1 + xi * scaled, -1 / xi) - ) - return bound(logc_expression, 1 + xi * scaled > 0, sigma > 0) - - def get_moment(value_var, size, mu, sigma, xi): - r""" - Using the mode, as the mean can be infinite when :math:`\xi > 1` - """ - mode = at.switch( - at.isclose(xi, 0), mu, mu + sigma * (at.pow(1 + xi, -xi) - 1) / xi - ) - return at.full(size, mode, dtype=aesara.config.floatX) diff --git a/pymc_experimental/tests/test_distributions.py b/pymc_experimental/tests/test_distributions.py index 479f9a18..78c7011e 100644 --- a/pymc_experimental/tests/test_distributions.py +++ b/pymc_experimental/tests/test_distributions.py @@ -11,51 +11,44 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -import functools -import itertools -import sys -import aesara -import aesara.tensor as at -import numpy as np -import numpy.random as nr +# general imports +import scipy.stats.distributions as ssd -import pytest -import scipy.stats -import scipy.stats.distributions as sp - -from aesara.compile.mode import Mode -from aesara.graph.basic import ancestors -from aesara.tensor.random.op import RandomVariable -from aesara.tensor.var import TensorVariable -from numpy import array, inf, log -from numpy.testing import assert_allclose, assert_almost_equal, assert_equal -from scipy import integrate -from scipy.special import erf, logit - -import pymc as pm +# test support imports from pymc +from pymc.test.test_distributions import ( + R, + Rplus, + Domain, + TestMatchesScipy, +) -from pymc.aesaraf import floatX, intX -from pymc-experimental.distributions import ( +# the distributions to be tested +from pymc_experimental.distributions import ( GenExtreme, ) -from pymc.distributions.shape_utils import to_tuple -from pymc.math import kronecker -from pymc.model import Deterministic, Model, Point, Potential -from pymc.tests.helpers import select_by_precision -from pymc.vartypes import continuous_types -def test_genextreme(self): - self.check_logp( - GenExtreme, - R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, - lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), - ) - self.check_logcdf( - GenExtreme, - R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, - lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), - ) - + +class TestMatchesScipyX(TestMatchesScipy): + """ + Wrapper class so that tests of experimental additions can be dropped into + PyMC directly on adoption. + """ + + def test_genextreme(self): + self.check_logp( + GenExtreme, + R, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + lambda value, mu, sigma, xi: ssd.genextreme.logpdf( + value, c=-xi, loc=mu, scale=sigma + ), + ) + self.check_logcdf( + GenExtreme, + R, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + lambda value, mu, sigma, xi: ssd.genextreme.logcdf( + value, c=-xi, loc=mu, scale=sigma + ), + ) diff --git a/pymc_experimental/tests/test_distributions_random.py b/pymc_experimental/tests/test_distributions_random.py index 76a5dfbb..540c2198 100644 --- a/pymc_experimental/tests/test_distributions_random.py +++ b/pymc_experimental/tests/test_distributions_random.py @@ -11,23 +11,21 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -import functools -import itertools -from typing import Callable, List, Optional +# general imports -import aesara -import numpy as np -import numpy.random as nr -import numpy.testing as npt -import pytest -import scipy.stats as st +# test support imports from pymc +from pymc.tests.test_distributions_random import ( + BaseTestDistribution, + seeded_scipy_distribution_builder, +) -from numpy.testing import assert_almost_equal, assert_array_almost_equal +# the distributions to be tested +import pymc_experimental as pmx class TestGenExtreme(BaseTestDistribution): - pymc_dist = pm.GenExtreme + pymc_dist = pmx.GenExtreme pymc_dist_params = {"mu": 0, "sigma": 1, "xi": -0.1} expected_rv_op_params = {"mu": 0, "sigma": 1, "xi": -0.1} # Notice, using different parametrization of xi sign to scipy From 8cfc35d49f9b34d2169ece09d572d792ad9014a7 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 6 Jul 2022 23:26:45 +1000 Subject: [PATCH 04/15] Added moment test --- pymc_experimental/tests/test_distributions.py | 10 ++-- .../tests/test_distributions_moments.py | 51 +++++++++++++++++++ 2 files changed, 54 insertions(+), 7 deletions(-) create mode 100644 pymc_experimental/tests/test_distributions_moments.py diff --git a/pymc_experimental/tests/test_distributions.py b/pymc_experimental/tests/test_distributions.py index 78c7011e..75ff123d 100644 --- a/pymc_experimental/tests/test_distributions.py +++ b/pymc_experimental/tests/test_distributions.py @@ -16,7 +16,7 @@ import scipy.stats.distributions as ssd # test support imports from pymc -from pymc.test.test_distributions import ( +from pymc.tests.test_distributions import ( R, Rplus, Domain, @@ -40,15 +40,11 @@ def test_genextreme(self): GenExtreme, R, {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, - lambda value, mu, sigma, xi: ssd.genextreme.logpdf( - value, c=-xi, loc=mu, scale=sigma - ), + lambda value, mu, sigma, xi: ssd.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), ) self.check_logcdf( GenExtreme, R, {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, - lambda value, mu, sigma, xi: ssd.genextreme.logcdf( - value, c=-xi, loc=mu, scale=sigma - ), + lambda value, mu, sigma, xi: ssd.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), ) diff --git a/pymc_experimental/tests/test_distributions_moments.py b/pymc_experimental/tests/test_distributions_moments.py new file mode 100644 index 00000000..483f24cf --- /dev/null +++ b/pymc_experimental/tests/test_distributions_moments.py @@ -0,0 +1,51 @@ +# import aesara +import numpy as np +import pytest +import scipy.stats as st + +# from aesara import tensor as at +# from scipy import special + +from pymc_experimental.distributions import ( + GenExtreme, +) + +from pymc.model import Model + +from pymc.tests.test_distributions_moments import assert_moment_is_expected + + +@pytest.mark.parametrize( + "mu, sigma, xi, size, expected", + [ + (0, 1, 0, None, 0), + (1, np.arange(1, 4), 0.1, None, np.arange(1, 4) * (1.1 ** -0.1 - 1) / 0.1), + (np.arange(5), 1, 0.1, None, np.arange(5) + (1.1 ** -0.1 - 1) / 0.1), + ( + 0, + 1, + np.linspace(-0.2, 0.2, 6), + None, + ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) + / np.linspace(-0.2, 0.2, 6), + ), + (1, 2, 0.1, 5, np.full(5, 1 + 2 * (1.1 ** -0.1 - 1) / 0.1)), + ( + np.arange(6), + np.arange(1, 7), + np.linspace(-0.2, 0.2, 6), + (3, 6), + np.full( + (3, 6), + np.arange(6) + + np.arange(1, 7) + * ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) + / np.linspace(-0.2, 0.2, 6), + ), + ), + ], +) +def test_genextreme_moment(mu, sigma, xi, size, expected): + with Model() as model: + GenExtreme("x", mu=mu, sigma=sigma, xi=xi, size=size) + assert_moment_is_expected(model, expected) From 85c3782744421e24f9bbfdb2f9b91f8fa466714d Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Tue, 27 Sep 2022 13:02:52 +1000 Subject: [PATCH 05/15] Update per #84 review; tests, and dist code --- pymc_experimental/distributions/continuous.py | 27 +++-- .../tests/distributions/__init__.py | 2 + .../tests/distributions/test_continuous.py | 108 ++++++++++++++++++ pymc_experimental/tests/test_distributions.py | 50 -------- .../tests/test_distributions_moments.py | 51 --------- .../tests/test_distributions_random.py | 38 ------ .../tests/test_histogram_approximation.py | 4 +- 7 files changed, 132 insertions(+), 148 deletions(-) create mode 100644 pymc_experimental/tests/distributions/__init__.py create mode 100644 pymc_experimental/tests/distributions/test_continuous.py delete mode 100644 pymc_experimental/tests/test_distributions.py delete mode 100644 pymc_experimental/tests/test_distributions_moments.py delete mode 100644 pymc_experimental/tests/test_distributions_random.py diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index 0b8e5246..d9129805 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -178,13 +178,17 @@ def logp(value, mu, sigma, xi): - ((xi + 1) / xi) * at.log1p(xi * scaled) - at.pow(1 + xi * scaled, -1 / xi), ) - # bnd = mu - sigma/xi - return check_parameters( - logp_expression, + # Confirm in valid domain + logp = at.switch( 1 + xi * (value - mu) / sigma > 0, - # at.switch(xi > 0, value > bnd, value < bnd), - sigma > 0, - ) + logp_expression, + -np.inf) + + return check_parameters( + logp, + sigma < 0, + 1 + xi * scaled < 0, + msg="sigma < 0, 1+xi*(x-mu)/sigma < 0") def logcdf(value, mu, sigma, xi): """ @@ -206,7 +210,16 @@ def logcdf(value, mu, sigma, xi): logc_expression = at.switch( at.isclose(xi, 0), -at.exp(-scaled), -at.pow(1 + xi * scaled, -1 / xi) ) - return check_parameters(logc_expression, 1 + xi * scaled > 0, sigma > 0) + # Confirm in valid domain + logc = at.switch( + 1 + xi * (value - mu) / sigma > 0, + logc_expression, + -np.inf) + + return check_parameters(logc, + sigma < 0, + 1 + xi * scaled < 0, + msg="sigma < 0, 1+xi*(x-mu)/sigma < 0") def get_moment(value_var, size, mu, sigma, xi): r""" diff --git a/pymc_experimental/tests/distributions/__init__.py b/pymc_experimental/tests/distributions/__init__.py new file mode 100644 index 00000000..41ff5483 --- /dev/null +++ b/pymc_experimental/tests/distributions/__init__.py @@ -0,0 +1,2 @@ +from pymc_experimental.distributions import histogram_utils +from pymc_experimental.distributions.histogram_utils import histogram_approximation diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py new file mode 100644 index 00000000..2373b34e --- /dev/null +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -0,0 +1,108 @@ + +# Copyright 2020 The PyMC Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# general imports +import numpy as np +import pytest +import scipy.stats as st +import scipy.stats.distributions as sp + +import pymc as pm + +# test support imports from pymc +from pymc.tests.distributions.util import ( + BaseTestDistributionRandom, + R, + Rplus, + Domain, + check_logp, + check_logcdf, + assert_moment_is_expected, + seeded_scipy_distribution_builder, +) + +# the distributions to be tested +from pymc_experimental.distributions import ( + GenExtreme, +) + +class TestGenExtreme(): + """ + Wrapper class so that tests of experimental additions can be dropped into + PyMC directly on adoption. + """ + + def test_genextreme(): + check_logp( + GenExtreme, + R, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), + ) + check_logcdf( + GenExtreme, + R, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), + ) + + @pytest.mark.parametrize( + "mu, sigma, xi, size, expected", + [ + (0, 1, 0, None, 0), + (1, np.arange(1, 4), 0.1, None, np.arange(1, 4) * (1.1 ** -0.1 - 1) / 0.1), + (np.arange(5), 1, 0.1, None, np.arange(5) + (1.1 ** -0.1 - 1) / 0.1), + ( + 0, + 1, + np.linspace(-0.2, 0.2, 6), + None, + ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) + / np.linspace(-0.2, 0.2, 6), + ), + (1, 2, 0.1, 5, np.full(5, 1 + 2 * (1.1 ** -0.1 - 1) / 0.1)), + ( + np.arange(6), + np.arange(1, 7), + np.linspace(-0.2, 0.2, 6), + (3, 6), + np.full( + (3, 6), + np.arange(6) + + np.arange(1, 7) + * ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) + / np.linspace(-0.2, 0.2, 6), + ), + ), + ], + ) + def test_genextreme_moment(mu, sigma, xi, size, expected): + with pm.Model() as model: + GenExtreme("x", mu=mu, sigma=sigma, xi=xi, size=size) + assert_moment_is_expected(model, expected) + + +class TestGenExtreme(BaseTestDistributionRandom): + pymc_dist = GenExtreme + pymc_dist_params = {"mu": 0, "sigma": 1, "xi": -0.1} + expected_rv_op_params = {"mu": 0, "sigma": 1, "xi": -0.1} + # Notice, using different parametrization of xi sign to scipy + reference_dist_params = {"loc": 0, "scale": 1, "c": 0.1} + reference_dist = seeded_scipy_distribution_builder("genextreme") + tests_to_run = [ + "check_pymc_params_match_rv_op", + "check_pymc_draws_match_reference", + "check_rv_size", + ] diff --git a/pymc_experimental/tests/test_distributions.py b/pymc_experimental/tests/test_distributions.py deleted file mode 100644 index 75ff123d..00000000 --- a/pymc_experimental/tests/test_distributions.py +++ /dev/null @@ -1,50 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -# general imports -import scipy.stats.distributions as ssd - -# test support imports from pymc -from pymc.tests.test_distributions import ( - R, - Rplus, - Domain, - TestMatchesScipy, -) - -# the distributions to be tested -from pymc_experimental.distributions import ( - GenExtreme, -) - - -class TestMatchesScipyX(TestMatchesScipy): - """ - Wrapper class so that tests of experimental additions can be dropped into - PyMC directly on adoption. - """ - - def test_genextreme(self): - self.check_logp( - GenExtreme, - R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, - lambda value, mu, sigma, xi: ssd.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), - ) - self.check_logcdf( - GenExtreme, - R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, - lambda value, mu, sigma, xi: ssd.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), - ) diff --git a/pymc_experimental/tests/test_distributions_moments.py b/pymc_experimental/tests/test_distributions_moments.py deleted file mode 100644 index 483f24cf..00000000 --- a/pymc_experimental/tests/test_distributions_moments.py +++ /dev/null @@ -1,51 +0,0 @@ -# import aesara -import numpy as np -import pytest -import scipy.stats as st - -# from aesara import tensor as at -# from scipy import special - -from pymc_experimental.distributions import ( - GenExtreme, -) - -from pymc.model import Model - -from pymc.tests.test_distributions_moments import assert_moment_is_expected - - -@pytest.mark.parametrize( - "mu, sigma, xi, size, expected", - [ - (0, 1, 0, None, 0), - (1, np.arange(1, 4), 0.1, None, np.arange(1, 4) * (1.1 ** -0.1 - 1) / 0.1), - (np.arange(5), 1, 0.1, None, np.arange(5) + (1.1 ** -0.1 - 1) / 0.1), - ( - 0, - 1, - np.linspace(-0.2, 0.2, 6), - None, - ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) - / np.linspace(-0.2, 0.2, 6), - ), - (1, 2, 0.1, 5, np.full(5, 1 + 2 * (1.1 ** -0.1 - 1) / 0.1)), - ( - np.arange(6), - np.arange(1, 7), - np.linspace(-0.2, 0.2, 6), - (3, 6), - np.full( - (3, 6), - np.arange(6) - + np.arange(1, 7) - * ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) - / np.linspace(-0.2, 0.2, 6), - ), - ), - ], -) -def test_genextreme_moment(mu, sigma, xi, size, expected): - with Model() as model: - GenExtreme("x", mu=mu, sigma=sigma, xi=xi, size=size) - assert_moment_is_expected(model, expected) diff --git a/pymc_experimental/tests/test_distributions_random.py b/pymc_experimental/tests/test_distributions_random.py deleted file mode 100644 index 540c2198..00000000 --- a/pymc_experimental/tests/test_distributions_random.py +++ /dev/null @@ -1,38 +0,0 @@ -# Copyright 2020 The PyMC Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -# general imports - -# test support imports from pymc -from pymc.tests.test_distributions_random import ( - BaseTestDistribution, - seeded_scipy_distribution_builder, -) - -# the distributions to be tested -import pymc_experimental as pmx - - -class TestGenExtreme(BaseTestDistribution): - pymc_dist = pmx.GenExtreme - pymc_dist_params = {"mu": 0, "sigma": 1, "xi": -0.1} - expected_rv_op_params = {"mu": 0, "sigma": 1, "xi": -0.1} - # Notice, using different parametrization of xi sign to scipy - reference_dist_params = {"loc": 0, "scale": 1, "c": 0.1} - reference_dist = seeded_scipy_distribution_builder("genextreme") - tests_to_run = [ - "check_pymc_params_match_rv_op", - "check_pymc_draws_match_reference", - "check_rv_size", - ] diff --git a/pymc_experimental/tests/test_histogram_approximation.py b/pymc_experimental/tests/test_histogram_approximation.py index e1eeda0e..33bda545 100644 --- a/pymc_experimental/tests/test_histogram_approximation.py +++ b/pymc_experimental/tests/test_histogram_approximation.py @@ -72,7 +72,7 @@ def test_histogram_approx_cont(use_dask, ndims): with pm.Model(): m = pm.Normal("m") s = pm.HalfNormal("s", size=2 if ndims > 1 else 1) - pot = pmx.distributions.histogram_approximation( + pot = pmx.distributions.histogram_utils.histogram_approximation( "histogram_potential", pm.Normal.dist(m, s), observed=data, n_quantiles=1000 ) trace = pm.sample(10, tune=0) # very fast @@ -88,7 +88,7 @@ def test_histogram_approx_discrete(use_dask, ndims): data = dask_df.from_array(data) with pm.Model(): s = pm.Exponential("s", 1.0, size=2 if ndims > 1 else 1) - pot = pmx.distributions.histogram_approximation( + pot = pmx.distributions.histogram_utils.histogram_approximation( "histogram_potential", pm.Poisson.dist(s), observed=data, min_count=10 ) trace = pm.sample(10, tune=0) # very fast From c43423f30e20632ed0d724e53ac4cadd85926dd6 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Tue, 27 Sep 2022 19:11:49 +1000 Subject: [PATCH 06/15] Updated tests and class per comments --- pymc_experimental/distributions/__init__.py | 2 +- pymc_experimental/distributions/continuous.py | 24 +++++++++---------- .../tests/distributions/test_continuous.py | 19 +++++++++++---- 3 files changed, 27 insertions(+), 18 deletions(-) diff --git a/pymc_experimental/distributions/__init__.py b/pymc_experimental/distributions/__init__.py index c04f53de..83302699 100644 --- a/pymc_experimental/distributions/__init__.py +++ b/pymc_experimental/distributions/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2020 The PyMC Developers +# Copyright 2022 The PyMC Developers # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index d9129805..20a0aad0 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -1,4 +1,4 @@ -# Copyright 2020 The PyMC Developers +# Copyright 2022 The PyMC Developers # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -19,7 +19,7 @@ The imports from pymc are not fully replicated here: add imports as necessary. """ -from typing import List, Tuple +from typing import List, Tuple, Union import aesara import aesara.tensor as at import numpy as np @@ -47,7 +47,7 @@ def __call__( @classmethod def rng_fn( cls, - rng: np.random.RandomState, + rng: Union[np.random.RandomState, np.random.Generator], mu: np.ndarray, sigma: np.ndarray, xi: np.ndarray, @@ -178,17 +178,17 @@ def logp(value, mu, sigma, xi): - ((xi + 1) / xi) * at.log1p(xi * scaled) - at.pow(1 + xi * scaled, -1 / xi), ) - # Confirm in valid domain + logp = at.switch( - 1 + xi * (value - mu) / sigma > 0, + at.gt((1 + xi * (value - mu) / sigma), 0.0), logp_expression, -np.inf) return check_parameters( logp, - sigma < 0, - 1 + xi * scaled < 0, - msg="sigma < 0, 1+xi*(x-mu)/sigma < 0") + sigma > 0, + 1 + xi * scaled > 0, + msg="sigma <= 0 or 1+xi*(x-mu)/sigma <= 0") def logcdf(value, mu, sigma, xi): """ @@ -210,16 +210,16 @@ def logcdf(value, mu, sigma, xi): logc_expression = at.switch( at.isclose(xi, 0), -at.exp(-scaled), -at.pow(1 + xi * scaled, -1 / xi) ) - # Confirm in valid domain + logc = at.switch( 1 + xi * (value - mu) / sigma > 0, logc_expression, -np.inf) return check_parameters(logc, - sigma < 0, - 1 + xi * scaled < 0, - msg="sigma < 0, 1+xi*(x-mu)/sigma < 0") + sigma > 0, + 1 + xi * scaled > 0, + msg="sigma <= 0 or 1+xi*(x-mu)/sigma <= 0") def get_moment(value_var, size, mu, sigma, xi): r""" diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index 2373b34e..939e64a1 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -38,23 +38,25 @@ GenExtreme, ) -class TestGenExtreme(): +class TestGenExtremeClass: """ Wrapper class so that tests of experimental additions can be dropped into PyMC directly on adoption. + + pm.logp(GenExtreme.dist(mu=0.,sigma=1.,xi=0.5),value=-0.01) """ - def test_genextreme(): + def test_genextreme(self): check_logp( GenExtreme, R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), ) check_logcdf( GenExtreme, R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -1, -0.5, 0, 0.5, 1, 1])}, + {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), ) @@ -88,11 +90,18 @@ def test_genextreme(): ), ], ) - def test_genextreme_moment(mu, sigma, xi, size, expected): + def test_genextreme_moment(self, mu, sigma, xi, size, expected): with pm.Model() as model: GenExtreme("x", mu=mu, sigma=sigma, xi=xi, size=size) assert_moment_is_expected(model, expected) + def test_gen_extreme_scipy_kwarg(self): + dist = GenExtreme.dist(xi=1, scipy=False) + assert dist.owner.inputs[-1].eval() == 1 + + dist = GenExtreme.dist(xi=1, scipy=True) + assert dist.owner.inputs[-1].eval() == -1 + class TestGenExtreme(BaseTestDistributionRandom): pymc_dist = GenExtreme From 13f9ade80a219fe1487035944b4a780b4ba9aa28 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 00:15:54 +1000 Subject: [PATCH 07/15] Tweaks post-review --- pymc_experimental/distributions/continuous.py | 2 +- pymc_experimental/tests/distributions/test_continuous.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index 20a0aad0..4ebb601b 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -180,7 +180,7 @@ def logp(value, mu, sigma, xi): ) logp = at.switch( - at.gt((1 + xi * (value - mu) / sigma), 0.0), + at.gt(1 + xi * scaled, 0.0), logp_expression, -np.inf) diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index 939e64a1..dcca994f 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -51,7 +51,8 @@ def test_genextreme(self): GenExtreme, R, {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, - lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma), + lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma) + if 1 + xi*(value-mu)/sigma > 0 else -np.inf ) check_logcdf( GenExtreme, From b181815a82e8cc89695da56376f4ff9f4384c637 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 03:10:08 +1000 Subject: [PATCH 08/15] Test tweaks and all pass --- pymc_experimental/distributions/continuous.py | 9 +++++++-- pymc_experimental/tests/distributions/test_continuous.py | 5 +++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index 4ebb601b..84476978 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -30,6 +30,7 @@ from aesara.tensor.var import TensorVariable from pymc.aesaraf import floatX from pymc.distributions.dist_math import check_parameters +from pymc.distributions.shape_utils import rv_size_is_none class GenExtremeRV(RandomVariable): @@ -188,6 +189,7 @@ def logp(value, mu, sigma, xi): logp, sigma > 0, 1 + xi * scaled > 0, + at.and_(xi > -1, xi < 1), msg="sigma <= 0 or 1+xi*(x-mu)/sigma <= 0") def logcdf(value, mu, sigma, xi): @@ -219,13 +221,16 @@ def logcdf(value, mu, sigma, xi): return check_parameters(logc, sigma > 0, 1 + xi * scaled > 0, + at.and_(xi > -1, xi < 1), msg="sigma <= 0 or 1+xi*(x-mu)/sigma <= 0") - def get_moment(value_var, size, mu, sigma, xi): + def moment(rv, size, mu, sigma, xi): r""" Using the mode, as the mean can be infinite when :math:`\xi > 1` """ mode = at.switch( at.isclose(xi, 0), mu, mu + sigma * (at.pow(1 + xi, -xi) - 1) / xi ) - return at.full(size, mode, dtype=aesara.config.floatX) + if not rv_size_is_none(size): + mode = at.full(size, mode) + return mode diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index dcca994f..13d7abb0 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -58,14 +58,15 @@ def test_genextreme(self): GenExtreme, R, {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, - lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma), + lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma) + if 1 + xi*(value-mu)/sigma > 0 else -np.inf ) @pytest.mark.parametrize( "mu, sigma, xi, size, expected", [ (0, 1, 0, None, 0), - (1, np.arange(1, 4), 0.1, None, np.arange(1, 4) * (1.1 ** -0.1 - 1) / 0.1), + (1, np.arange(1, 4), 0.1, None, 1 + np.arange(1, 4) * (1.1 ** -0.1 - 1) / 0.1), (np.arange(5), 1, 0.1, None, np.arange(5) + (1.1 ** -0.1 - 1) / 0.1), ( 0, From 14c825202dafcae8ae78747fc76d783a6929d5bc Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 03:51:24 +1000 Subject: [PATCH 09/15] Resolve unnecessary check --- pymc_experimental/distributions/continuous.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index 84476978..208e18c5 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -188,9 +188,8 @@ def logp(value, mu, sigma, xi): return check_parameters( logp, sigma > 0, - 1 + xi * scaled > 0, at.and_(xi > -1, xi < 1), - msg="sigma <= 0 or 1+xi*(x-mu)/sigma <= 0") + msg="sigma <= 0 or -1 < xi < 1") def logcdf(value, mu, sigma, xi): """ @@ -220,9 +219,8 @@ def logcdf(value, mu, sigma, xi): return check_parameters(logc, sigma > 0, - 1 + xi * scaled > 0, at.and_(xi > -1, xi < 1), - msg="sigma <= 0 or 1+xi*(x-mu)/sigma <= 0") + msg="sigma <= 0 or -1 < xi < 1") def moment(rv, size, mu, sigma, xi): r""" From 8471fcad9d6825e30e041038c90033570c1a9bc0 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 11:38:13 +1000 Subject: [PATCH 10/15] Reverse constraints msg --- pymc_experimental/distributions/continuous.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index 208e18c5..33abadfe 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -189,7 +189,7 @@ def logp(value, mu, sigma, xi): logp, sigma > 0, at.and_(xi > -1, xi < 1), - msg="sigma <= 0 or -1 < xi < 1") + msg="sigma > 0 or -1 < xi < 1") def logcdf(value, mu, sigma, xi): """ @@ -220,7 +220,7 @@ def logcdf(value, mu, sigma, xi): return check_parameters(logc, sigma > 0, at.and_(xi > -1, xi < 1), - msg="sigma <= 0 or -1 < xi < 1") + msg="sigma > 0 or -1 < xi < 1") def moment(rv, size, mu, sigma, xi): r""" From e90be67dd9251e6e111a5d8bc143cd7f2a3fcbcf Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 12:11:29 +1000 Subject: [PATCH 11/15] Clearing precommit requirements --- pymc_experimental/distributions/__init__.py | 4 +- pymc_experimental/distributions/continuous.py | 42 ++++++------------- .../tests/distributions/test_continuous.py | 28 ++++++------- 3 files changed, 27 insertions(+), 47 deletions(-) diff --git a/pymc_experimental/distributions/__init__.py b/pymc_experimental/distributions/__init__.py index 83302699..fe7c7037 100644 --- a/pymc_experimental/distributions/__init__.py +++ b/pymc_experimental/distributions/__init__.py @@ -17,9 +17,7 @@ Experimental probability distributions for stochastic nodes in PyMC. """ -from pymc_experimental.distributions.continuous import ( - GenExtreme, -) +from pymc_experimental.distributions.continuous import GenExtreme __all__ = [ "GenExtreme", diff --git a/pymc_experimental/distributions/continuous.py b/pymc_experimental/distributions/continuous.py index 33abadfe..185cb94a 100644 --- a/pymc_experimental/distributions/continuous.py +++ b/pymc_experimental/distributions/continuous.py @@ -20,17 +20,16 @@ """ from typing import List, Tuple, Union -import aesara + import aesara.tensor as at import numpy as np -from scipy import stats - from aesara.tensor.random.op import RandomVariable -from pymc.distributions.distribution import Continuous from aesara.tensor.var import TensorVariable from pymc.aesaraf import floatX from pymc.distributions.dist_math import check_parameters +from pymc.distributions.distribution import Continuous from pymc.distributions.shape_utils import rv_size_is_none +from scipy import stats class GenExtremeRV(RandomVariable): @@ -40,9 +39,7 @@ class GenExtremeRV(RandomVariable): dtype: str = "floatX" _print_name: Tuple[str, str] = ("Generalized Extreme Value", "\\operatorname{GEV}") - def __call__( - self, mu=0.0, sigma=1.0, xi=0.0, size=None, **kwargs - ) -> TensorVariable: + def __call__(self, mu=0.0, sigma=1.0, xi=0.0, size=None, **kwargs) -> TensorVariable: return super().__call__(mu, sigma, xi, size=size, **kwargs) @classmethod @@ -55,9 +52,7 @@ def rng_fn( size: Tuple[int, ...], ) -> np.ndarray: # Notice negative here, since remainder of GenExtreme is based on Coles parametrization - return stats.genextreme.rvs( - c=-xi, loc=mu, scale=sigma, random_state=rng, size=size - ) + return stats.genextreme.rvs(c=-xi, loc=mu, scale=sigma, random_state=rng, size=size) gev = GenExtremeRV() @@ -180,16 +175,11 @@ def logp(value, mu, sigma, xi): - at.pow(1 + xi * scaled, -1 / xi), ) - logp = at.switch( - at.gt(1 + xi * scaled, 0.0), - logp_expression, - -np.inf) + logp = at.switch(at.gt(1 + xi * scaled, 0.0), logp_expression, -np.inf) return check_parameters( - logp, - sigma > 0, - at.and_(xi > -1, xi < 1), - msg="sigma > 0 or -1 < xi < 1") + logp, sigma > 0, at.and_(xi > -1, xi < 1), msg="sigma > 0 or -1 < xi < 1" + ) def logcdf(value, mu, sigma, xi): """ @@ -212,23 +202,17 @@ def logcdf(value, mu, sigma, xi): at.isclose(xi, 0), -at.exp(-scaled), -at.pow(1 + xi * scaled, -1 / xi) ) - logc = at.switch( - 1 + xi * (value - mu) / sigma > 0, - logc_expression, - -np.inf) + logc = at.switch(1 + xi * (value - mu) / sigma > 0, logc_expression, -np.inf) - return check_parameters(logc, - sigma > 0, - at.and_(xi > -1, xi < 1), - msg="sigma > 0 or -1 < xi < 1") + return check_parameters( + logc, sigma > 0, at.and_(xi > -1, xi < 1), msg="sigma > 0 or -1 < xi < 1" + ) def moment(rv, size, mu, sigma, xi): r""" Using the mode, as the mean can be infinite when :math:`\xi > 1` """ - mode = at.switch( - at.isclose(xi, 0), mu, mu + sigma * (at.pow(1 + xi, -xi) - 1) / xi - ) + mode = at.switch(at.isclose(xi, 0), mu, mu + sigma * (at.pow(1 + xi, -xi) - 1) / xi) if not rv_size_is_none(size): mode = at.full(size, mode) return mode diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index 13d7abb0..80fd6ae3 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -1,4 +1,3 @@ - # Copyright 2020 The PyMC Developers # # Licensed under the Apache License, Version 2.0 (the "License"); @@ -15,28 +14,25 @@ # general imports import numpy as np +import pymc as pm import pytest -import scipy.stats as st import scipy.stats.distributions as sp -import pymc as pm - # test support imports from pymc from pymc.tests.distributions.util import ( BaseTestDistributionRandom, + Domain, R, Rplus, - Domain, - check_logp, - check_logcdf, assert_moment_is_expected, + check_logcdf, + check_logp, seeded_scipy_distribution_builder, ) # the distributions to be tested -from pymc_experimental.distributions import ( - GenExtreme, -) +from pymc_experimental.distributions import GenExtreme + class TestGenExtremeClass: """ @@ -52,22 +48,24 @@ def test_genextreme(self): R, {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma) - if 1 + xi*(value-mu)/sigma > 0 else -np.inf + if 1 + xi * (value - mu) / sigma > 0 + else -np.inf, ) check_logcdf( GenExtreme, R, {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma) - if 1 + xi*(value-mu)/sigma > 0 else -np.inf + if 1 + xi * (value - mu) / sigma > 0 + else -np.inf, ) @pytest.mark.parametrize( "mu, sigma, xi, size, expected", [ (0, 1, 0, None, 0), - (1, np.arange(1, 4), 0.1, None, 1 + np.arange(1, 4) * (1.1 ** -0.1 - 1) / 0.1), - (np.arange(5), 1, 0.1, None, np.arange(5) + (1.1 ** -0.1 - 1) / 0.1), + (1, np.arange(1, 4), 0.1, None, 1 + np.arange(1, 4) * (1.1**-0.1 - 1) / 0.1), + (np.arange(5), 1, 0.1, None, np.arange(5) + (1.1**-0.1 - 1) / 0.1), ( 0, 1, @@ -76,7 +74,7 @@ def test_genextreme(self): ((1 + np.linspace(-0.2, 0.2, 6)) ** -np.linspace(-0.2, 0.2, 6) - 1) / np.linspace(-0.2, 0.2, 6), ), - (1, 2, 0.1, 5, np.full(5, 1 + 2 * (1.1 ** -0.1 - 1) / 0.1)), + (1, 2, 0.1, 5, np.full(5, 1 + 2 * (1.1**-0.1 - 1) / 0.1)), ( np.arange(6), np.arange(1, 7), From 306ce90d0f94fa6452ce5eba6fb5e617e5e9409e Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 16:58:25 +1000 Subject: [PATCH 12/15] Skip tests on float32 due to underflow of pymc vs scipy --- pymc_experimental/tests/distributions/test_continuous.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index 80fd6ae3..4491186f 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -13,6 +13,7 @@ # limitations under the License. # general imports +import aesara import numpy as np import pymc as pm import pytest @@ -42,6 +43,10 @@ class TestGenExtremeClass: pm.logp(GenExtreme.dist(mu=0.,sigma=1.,xi=0.5),value=-0.01) """ + @pytest.mark.skipif( + condition=(aesara.config.floatX == "float32"), + reason="PyMC underflows earlier than scipy on float32", + ) def test_genextreme(self): check_logp( GenExtreme, @@ -51,6 +56,7 @@ def test_genextreme(self): if 1 + xi * (value - mu) / sigma > 0 else -np.inf, ) + check_logcdf( GenExtreme, R, @@ -59,6 +65,8 @@ def test_genextreme(self): if 1 + xi * (value - mu) / sigma > 0 else -np.inf, ) + if aesara.config.floatX == "float32": + raise Exception("Flaky test: It passed this time, but XPASS is not allowed.") @pytest.mark.parametrize( "mu, sigma, xi, size, expected", From 23010155348f97c72b156f357f67529b06fe0101 Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Wed, 28 Sep 2022 23:38:35 +1000 Subject: [PATCH 13/15] Improve test grid for parameters --- .../tests/distributions/test_continuous.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index 4491186f..bb73f47b 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -24,7 +24,7 @@ BaseTestDistributionRandom, Domain, R, - Rplus, + Rplusbig, assert_moment_is_expected, check_logcdf, check_logp, @@ -51,7 +51,11 @@ def test_genextreme(self): check_logp( GenExtreme, R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, + { + "mu": Domain([-np.inf, -2.1, -1, -0.01, 0.0, 0.01, 1, 2.1, np.inf]), + "sigma": Rplusbig, + "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1]), + }, lambda value, mu, sigma, xi: sp.genextreme.logpdf(value, c=-xi, loc=mu, scale=sigma) if 1 + xi * (value - mu) / sigma > 0 else -np.inf, @@ -60,7 +64,11 @@ def test_genextreme(self): check_logcdf( GenExtreme, R, - {"mu": R, "sigma": Rplus, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1])}, + { + "mu": Domain([-np.inf, -2.1, -1, -0.01, 0.0, 0.01, 1, 2.1, np.inf]), + "sigma": Rplusbig, + "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1]), + }, lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma) if 1 + xi * (value - mu) / sigma > 0 else -np.inf, From 0022a2f843d49f9bccae2a812adc9cf09a95bd6d Mon Sep 17 00:00:00 2001 From: Colin Caprani Date: Thu, 29 Sep 2022 00:05:34 +1000 Subject: [PATCH 14/15] Revert mu test domain to R --- pymc_experimental/tests/distributions/test_continuous.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index bb73f47b..add13743 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -52,7 +52,7 @@ def test_genextreme(self): GenExtreme, R, { - "mu": Domain([-np.inf, -2.1, -1, -0.01, 0.0, 0.01, 1, 2.1, np.inf]), + "mu": R, "sigma": Rplusbig, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1]), }, @@ -65,7 +65,7 @@ def test_genextreme(self): GenExtreme, R, { - "mu": Domain([-np.inf, -2.1, -1, -0.01, 0.0, 0.01, 1, 2.1, np.inf]), + "mu": R, "sigma": Rplusbig, "xi": Domain([-1, -0.99, -0.5, 0, 0.5, 0.99, 1]), }, From 544cf423690178f79f745d31e818d952a74d933c Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Wed, 28 Sep 2022 17:15:17 +0200 Subject: [PATCH 15/15] Separate logp and logcdf tests --- .../tests/distributions/test_continuous.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/pymc_experimental/tests/distributions/test_continuous.py b/pymc_experimental/tests/distributions/test_continuous.py index add13743..f87eaf0e 100644 --- a/pymc_experimental/tests/distributions/test_continuous.py +++ b/pymc_experimental/tests/distributions/test_continuous.py @@ -30,6 +30,7 @@ check_logp, seeded_scipy_distribution_builder, ) +from pymc.tests.helpers import select_by_precision # the distributions to be tested from pymc_experimental.distributions import GenExtreme @@ -43,11 +44,11 @@ class TestGenExtremeClass: pm.logp(GenExtreme.dist(mu=0.,sigma=1.,xi=0.5),value=-0.01) """ - @pytest.mark.skipif( + @pytest.mark.xfail( condition=(aesara.config.floatX == "float32"), reason="PyMC underflows earlier than scipy on float32", ) - def test_genextreme(self): + def test_logp(self): check_logp( GenExtreme, R, @@ -61,6 +62,10 @@ def test_genextreme(self): else -np.inf, ) + if aesara.config.floatX == "float32": + raise Exception("Flaky test: It passed this time, but XPASS is not allowed.") + + def test_logcdf(self): check_logcdf( GenExtreme, R, @@ -72,9 +77,8 @@ def test_genextreme(self): lambda value, mu, sigma, xi: sp.genextreme.logcdf(value, c=-xi, loc=mu, scale=sigma) if 1 + xi * (value - mu) / sigma > 0 else -np.inf, + decimal=select_by_precision(float64=6, float32=2), ) - if aesara.config.floatX == "float32": - raise Exception("Flaky test: It passed this time, but XPASS is not allowed.") @pytest.mark.parametrize( "mu, sigma, xi, size, expected",