From 613857dc687fd1bae592021c086ad111675d4f06 Mon Sep 17 00:00:00 2001 From: Trey Wenger Date: Wed, 27 Sep 2023 10:06:26 -0500 Subject: [PATCH 1/3] Update PyTensor dependency --- conda-envs/environment-dev.yml | 2 +- conda-envs/environment-docs.yml | 2 +- conda-envs/environment-test.yml | 2 +- conda-envs/windows-environment-dev.yml | 2 +- conda-envs/windows-environment-test.yml | 2 +- pymc/distributions/continuous.py | 25 +++++---- requirements-dev.txt | 2 +- requirements.txt | 2 +- tests/distributions/test_continuous.py | 2 +- tests/logprob/test_mixture.py | 72 +++++++++++-------------- 10 files changed, 52 insertions(+), 61 deletions(-) diff --git a/conda-envs/environment-dev.yml b/conda-envs/environment-dev.yml index db607bcfb6..b9a533a464 100644 --- a/conda-envs/environment-dev.yml +++ b/conda-envs/environment-dev.yml @@ -14,7 +14,7 @@ dependencies: - numpy>=1.15.0 - pandas>=0.24.0 - pip -- pytensor>=2.16.1,<2.17 +- pytensor>=2.17.0,<2.18 - python-graphviz - networkx - scipy>=1.4.1 diff --git a/conda-envs/environment-docs.yml b/conda-envs/environment-docs.yml index e42e20307b..e35839fbe8 100644 --- a/conda-envs/environment-docs.yml +++ b/conda-envs/environment-docs.yml @@ -12,7 +12,7 @@ dependencies: - numpy>=1.15.0 - pandas>=0.24.0 - pip -- pytensor>=2.16.1,<2.17 +- pytensor>=2.17.0,<2.18 - python-graphviz - scipy>=1.4.1 - typing-extensions>=3.7.4 diff --git a/conda-envs/environment-test.yml b/conda-envs/environment-test.yml index 1f26a72b11..01ad9c4f31 100644 --- a/conda-envs/environment-test.yml +++ b/conda-envs/environment-test.yml @@ -17,7 +17,7 @@ dependencies: - numpy>=1.15.0 - pandas>=0.24.0 - pip -- pytensor>=2.16.1,<2.17 +- pytensor>=2.17.0,<2.18 - python-graphviz - networkx - scipy>=1.4.1 diff --git a/conda-envs/windows-environment-dev.yml b/conda-envs/windows-environment-dev.yml index 40f0cb6c53..2c86d9af49 100644 --- a/conda-envs/windows-environment-dev.yml +++ b/conda-envs/windows-environment-dev.yml @@ -14,7 +14,7 @@ dependencies: - numpy>=1.15.0 - pandas>=0.24.0 - pip -- pytensor>=2.16.1,<2.17 +- pytensor>=2.17.0,<2.18 - python-graphviz - networkx - scipy>=1.4.1 diff --git a/conda-envs/windows-environment-test.yml b/conda-envs/windows-environment-test.yml index 9471a0d11b..9a0699752f 100644 --- a/conda-envs/windows-environment-test.yml +++ b/conda-envs/windows-environment-test.yml @@ -17,7 +17,7 @@ dependencies: - numpy>=1.15.0 - pandas>=0.24.0 - pip -- pytensor>=2.16.1,<2.17 +- pytensor>=2.17.0,<2.18 - python-graphviz - networkx - scipy>=1.4.1 diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 1b096667ae..4212e4e4dc 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -36,10 +36,10 @@ from pytensor.tensor.math import tanh from pytensor.tensor.random.basic import ( BetaRV, + _gamma, cauchy, chisquare, exponential, - gamma, gumbel, halfcauchy, halfnormal, @@ -2201,16 +2201,17 @@ class Gamma(PositiveContinuous): sigma : tensor_like of float, optional Alternative scale parameter (sigma > 0). """ - rv_op = gamma + # gamma is temporarily a deprecation wrapper in PyTensor + rv_op = _gamma @classmethod def dist(cls, alpha=None, beta=None, mu=None, sigma=None, **kwargs): alpha, beta = cls.get_alpha_beta(alpha, beta, mu, sigma) alpha = pt.as_tensor_variable(floatX(alpha)) beta = pt.as_tensor_variable(floatX(beta)) - - # The PyTensor `GammaRV` `Op` will invert the `beta` parameter itself - return super().dist([alpha, beta], **kwargs) + # PyTensor gamma op is parametrized in terms of scale (1/beta) + scale = pt.reciprocal(beta) + return super().dist([alpha, scale], **kwargs) @classmethod def get_alpha_beta(cls, alpha=None, beta=None, mu=None, sigma=None): @@ -2232,15 +2233,14 @@ def get_alpha_beta(cls, alpha=None, beta=None, mu=None, sigma=None): return alpha, beta - def moment(rv, size, alpha, inv_beta): - # The PyTensor `GammaRV` `Op` inverts the `beta` parameter itself - mean = alpha * inv_beta + def moment(rv, size, alpha, scale): + mean = alpha * scale if not rv_size_is_none(size): mean = pt.full(size, mean) return mean - def logp(value, alpha, inv_beta): - beta = pt.reciprocal(inv_beta) + def logp(value, alpha, scale): + beta = pt.reciprocal(scale) res = -pt.gammaln(alpha) + logpow(beta, alpha) - beta * value + logpow(value, alpha - 1) res = pt.switch(pt.ge(value, 0.0), res, -np.inf) return check_parameters( @@ -2250,14 +2250,13 @@ def logp(value, alpha, inv_beta): msg="alpha > 0, beta > 0", ) - def logcdf(value, alpha, inv_beta): - beta = pt.reciprocal(inv_beta) + def logcdf(value, alpha, scale): + beta = pt.reciprocal(scale) res = pt.switch( pt.lt(value, 0), -np.inf, pt.log(pt.gammainc(alpha, beta * value)), ) - return check_parameters(res, 0 < alpha, 0 < beta, msg="alpha > 0, beta > 0") diff --git a/requirements-dev.txt b/requirements-dev.txt index 71aff41b3b..18eede6a60 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -18,7 +18,7 @@ numpydoc pandas>=0.24.0 polyagamma pre-commit>=2.8.0 -pytensor>=2.16.1,<2.17 +pytensor>=2.17.0,<2.18 pytest-cov>=2.5 pytest>=3.0 scipy>=1.4.1 diff --git a/requirements.txt b/requirements.txt index 0cd2d7201d..183908fe5e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,6 +4,6 @@ cloudpickle fastprogress>=0.2.0 numpy>=1.15.0 pandas>=0.24.0 -pytensor>=2.16.1,<2.17 +pytensor>=2.17.0,<2.18 scipy>=1.4.1 typing-extensions>=3.7.4 diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index bab4e281c1..3d74b00bff 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -2199,7 +2199,7 @@ class TestHalfCauchy(BaseTestDistributionRandom): class TestGamma(BaseTestDistributionRandom): pymc_dist = pm.Gamma pymc_dist_params = {"alpha": 2.0, "beta": 5.0} - expected_rv_op_params = {"alpha": 2.0, "beta": 1 / 5.0} + expected_rv_op_params = {"shape": 2.0, "scale": 1 / 5.0} reference_dist_params = {"shape": 2.0, "scale": 1 / 5.0} reference_dist = seeded_numpy_distribution_builder("gamma") checks_to_run = [ diff --git a/tests/logprob/test_mixture.py b/tests/logprob/test_mixture.py index ef46b2cf27..1d2801592a 100644 --- a/tests/logprob/test_mixture.py +++ b/tests/logprob/test_mixture.py @@ -62,16 +62,14 @@ def test_mixture_basics(): - srng = pt.random.RandomStream(29833) - def create_mix_model(size, axis): - X_rv = srng.normal(0, 1, size=size, name="X") - Y_rv = srng.gamma(0.5, 0.5, size=size, name="Y") + X_rv = pt.random.normal(0, 1, size=size, name="X") + Y_rv = pt.random.gamma(0.5, scale=2.0, size=size, name="Y") p_at = pt.scalar("p") p_at.tag.test_value = 0.5 - I_rv = srng.bernoulli(p_at, size=size, name="I") + I_rv = pt.random.bernoulli(p_at, size=size, name="I") i_vv = I_rv.clone() i_vv.name = "i" @@ -119,15 +117,13 @@ def create_mix_model(size, axis): ], ) def test_compute_test_value(op_constructor): - srng = pt.random.RandomStream(29833) - - X_rv = srng.normal(0, 1, name="X") - Y_rv = srng.gamma(0.5, 0.5, name="Y") + X_rv = pt.random.normal(0, 1, name="X") + Y_rv = pt.random.gamma(0.5, scale=2.0, name="Y") p_at = pt.scalar("p") p_at.tag.test_value = 0.3 - I_rv = srng.bernoulli(p_at, name="I") + I_rv = pt.random.bernoulli(p_at, name="I") i_vv = I_rv.clone() i_vv.name = "i" @@ -160,20 +156,18 @@ def test_compute_test_value(op_constructor): ], ) def test_hetero_mixture_binomial(p_val, size, supported): - srng = pt.random.RandomStream(29833) - - X_rv = srng.normal(0, 1, size=size, name="X") - Y_rv = srng.gamma(0.5, 0.5, size=size, name="Y") + X_rv = pt.random.normal(0, 1, size=size, name="X") + Y_rv = pt.random.gamma(0.5, scale=2.0, size=size, name="Y") if np.ndim(p_val) == 0: p_at = pt.scalar("p") p_at.tag.test_value = p_val - I_rv = srng.bernoulli(p_at, size=size, name="I") + I_rv = pt.random.bernoulli(p_at, size=size, name="I") p_val_1 = p_val else: p_at = pt.vector("p") p_at.tag.test_value = np.array(p_val, dtype=pytensor.config.floatX) - I_rv = srng.categorical(p_at, size=size, name="I") + I_rv = pt.random.categorical(p_at, size=size, name="I") p_val_1 = p_val[1] i_vv = I_rv.clone() @@ -203,7 +197,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): bern_sp = sp.bernoulli(p_val_1) norm_sp = sp.norm(loc=0, scale=1) - gamma_sp = sp.gamma(0.5, scale=1.0 / 0.5) + gamma_sp = sp.gamma(0.5, scale=2.0) for i in range(10): i_val = bern_sp.rvs(size=size, random_state=test_val_rng) @@ -230,7 +224,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -251,7 +245,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array([0.5], dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array([100], dtype=pytensor.config.floatX), @@ -272,7 +266,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array([0.5], dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array([100], dtype=pytensor.config.floatX), @@ -293,7 +287,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -314,7 +308,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -335,7 +329,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -361,7 +355,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -384,7 +378,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -405,7 +399,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -426,7 +420,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -447,7 +441,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array([0.5], dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array([100], dtype=pytensor.config.floatX), @@ -468,7 +462,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array([0.5, 1], dtype=pytensor.config.floatX), - np.array([0.5, 1], dtype=pytensor.config.floatX), + np.array([2.0, 1], dtype=pytensor.config.floatX), ), ( np.array([100, 1000], dtype=pytensor.config.floatX), @@ -489,7 +483,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array([0.5, 1], dtype=pytensor.config.floatX), - np.array([0.5, 1], dtype=pytensor.config.floatX), + np.array([2.0, 1], dtype=pytensor.config.floatX), ), ( np.array([100, 1000], dtype=pytensor.config.floatX), @@ -510,7 +504,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -531,7 +525,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -552,7 +546,7 @@ def test_hetero_mixture_binomial(p_val, size, supported): ), ( np.array(0.5, dtype=pytensor.config.floatX), - np.array(0.5, dtype=pytensor.config.floatX), + np.array(2.0, dtype=pytensor.config.floatX), ), ( np.array(100, dtype=pytensor.config.floatX), @@ -570,16 +564,14 @@ def test_hetero_mixture_binomial(p_val, size, supported): def test_hetero_mixture_categorical( X_args, Y_args, Z_args, p_val, comp_size, idx_size, extra_indices, join_axis, supported ): - srng = pt.random.RandomStream(29833) - - X_rv = srng.normal(*X_args, size=comp_size, name="X") - Y_rv = srng.gamma(*Y_args, size=comp_size, name="Y") - Z_rv = srng.normal(*Z_args, size=comp_size, name="Z") + X_rv = pt.random.normal(*X_args, size=comp_size, name="X") + Y_rv = pt.random.gamma(Y_args[0], scale=Y_args[1], size=comp_size, name="Y") + Z_rv = pt.random.normal(*Z_args, size=comp_size, name="Z") p_at = pt.as_tensor(p_val).type() p_at.name = "p" p_at.tag.test_value = np.array(p_val, dtype=pytensor.config.floatX) - I_rv = srng.categorical(p_at, size=idx_size, name="I") + I_rv = pt.random.categorical(p_at, size=idx_size, name="I") i_vv = I_rv.clone() i_vv.name = "i" @@ -612,7 +604,7 @@ def test_hetero_mixture_categorical( test_val_rng = np.random.RandomState(3238) norm_1_sp = sp.norm(loc=X_args[0], scale=X_args[1]) - gamma_sp = sp.gamma(Y_args[0], scale=1 / Y_args[1]) + gamma_sp = sp.gamma(Y_args[0], scale=Y_args[1]) norm_2_sp = sp.norm(loc=Z_args[0], scale=Z_args[1]) # Handle scipy annoying squeeze of random draws From c6fca17f94701ee2eab61f5851b40760353cc543 Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Mon, 2 Oct 2023 12:37:32 +0100 Subject: [PATCH 2/3] Add regression test for Truncated Gamma --- tests/distributions/test_truncated.py | 32 ++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/tests/distributions/test_truncated.py b/tests/distributions/test_truncated.py index 53e933285b..f4ddfef1e9 100644 --- a/tests/distributions/test_truncated.py +++ b/tests/distributions/test_truncated.py @@ -20,7 +20,7 @@ from pytensor.tensor.random.basic import GeometricRV, NormalRV from pymc import Censored, Model, draw, find_MAP -from pymc.distributions.continuous import Exponential, TruncatedNormalRV +from pymc.distributions.continuous import Exponential, Gamma, TruncatedNormalRV from pymc.distributions.shape_utils import change_dist_size from pymc.distributions.transforms import _default_transform from pymc.distributions.truncated import Truncated, TruncatedRV, _truncated @@ -392,3 +392,33 @@ def test_truncated_inference(): map = find_MAP(progressbar=False) assert np.isclose(map["lam"], lam_true, atol=0.1) + + +def test_truncated_gamma(): + # Regression test for https://github.com/pymc-devs/pymc/issues/6931 + alpha = 3.0 + beta = 3.0 + upper = 2.5 + x = np.linspace(0.0, upper + 0.5, 100) + + gamma_scipy = scipy.stats.gamma(a=alpha, scale=1.0 / beta) + logp_scipy = gamma_scipy.logpdf(x) - gamma_scipy.logcdf(upper) + logp_scipy[x > upper] = -np.inf + + gamma_trunc_pymc = Truncated.dist( + Gamma.dist(alpha=alpha, beta=beta), + upper=upper, + ) + logp_pymc = logp(gamma_trunc_pymc, x).eval() + np.testing.assert_allclose( + logp_pymc, + logp_scipy, + ) + + # Changing the size used to invert the beta Gamma parameter again + resized_gamma_trunc_pymc = change_dist_size(gamma_trunc_pymc, new_size=x.shape) + logp_resized_pymc = logp(resized_gamma_trunc_pymc, x).eval() + np.testing.assert_allclose( + logp_resized_pymc, + logp_scipy, + ) From 1f8e470b7e140c76c3b3ef9456587de4187f371a Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Mon, 2 Oct 2023 12:49:25 +0100 Subject: [PATCH 3/3] Fix mypy failure --- pymc/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc/data.py b/pymc/data.py index d75fd698b8..0b18f2232a 100644 --- a/pymc/data.py +++ b/pymc/data.py @@ -234,7 +234,7 @@ def determine_coords( for dim in dims: dim_name = dim # str is applied because dim entries may be None - coords[str(dim_name)] = value[dim].to_numpy() + coords[str(dim_name)] = cast(xr.DataArray, value[dim]).to_numpy() if isinstance(value, np.ndarray) and dims is not None: if len(dims) != value.ndim: