From ecb755a1845e804ea143073ec9bae0df0995d90a Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Fri, 17 Feb 2023 11:42:12 +0100 Subject: [PATCH 1/9] added 3 new icdf functions and their corresponding tests --- pymc/distributions/continuous.py | 9 ++++++ pymc/tests/distributions/test_continuous.py | 31 +++++++++++++++++++-- 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index c8c874267f..0e3baecc0f 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -343,6 +343,9 @@ def logcdf(value, lower, upper): msg="lower <= upper", ) + def icdf(value, lower, upper): + return lower + (upper - lower) * value + @_default_transform.register(Uniform) def uniform_default_transform(op, rv): @@ -835,6 +838,9 @@ def logcdf(value, loc, sigma): msg="sigma > 0", ) + def icdf(value, lower, upper, loc, sigma): + return stats.truncnorm.ppf(q=value, a=lower, b=upper, loc=loc, scale=sigma) + class WaldRV(RandomVariable): name = "wald" @@ -1013,6 +1019,9 @@ def logcdf(value, mu, lam, alpha): msg="mu > 0, lam > 0, alpha >= 0", ) + def icdf(mu, lam): + return stats.invgauss(mu=mu, scale=lam) + class BetaClippedRV(BetaRV): @classmethod diff --git a/pymc/tests/distributions/test_continuous.py b/pymc/tests/distributions/test_continuous.py index 61cca44858..6745a65e75 100644 --- a/pymc/tests/distributions/test_continuous.py +++ b/pymc/tests/distributions/test_continuous.py @@ -26,7 +26,14 @@ import pymc as pm -from pymc.distributions.continuous import Normal, get_tau_sigma, interpolated +from pymc.distributions.continuous import ( + Normal, + TruncatedNormal, + Uniform, + Wald, + get_tau_sigma, + interpolated, +) from pymc.distributions.dist_math import clipped_beta_rvs from pymc.logprob.abstract import logcdf from pymc.logprob.joint_logprob import logp @@ -2281,5 +2288,25 @@ def test_normal_icdf(self, dist_params, obs, size): dist_params = dict(zip(dist_params_at, dist_params)) x = Normal.dist(*dist_params_at, size=size_at) - scipy_logprob_tester(x, obs, dist_params, test_fn=st.norm.ppf, test="icdf") + + def test_uniform_icdf(self, dist_params, obs, size): + dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) + dist_params = dict(zip(dist_params_at, dist_params)) + + x = Uniform.dist(*dist_params_at, size=size_at) + scipy_logprob_tester(x, obs, dist_params, test_fn=st.uniform.ppf, test="icdf") + + def test_truncated_normal_icdf(self, dist_params, obs, size): + dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) + dist_params = dict(zip(dist_params_at, dist_params)) + + x = TruncatedNormal.dist(*dist_params_at, size=size_at) + scipy_logprob_tester(x, obs, dist_params, test_fn=st.truncnorm.ppf, test="icdf") + + def test_wald_icdf(self, *dist_params, obs, size): + dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) + dist_params = dict(zip(dist_params_at, dist_params)) + + x = Wald.dist(*dist_params_at, size=size_at) + scipy_logprob_tester(x, obs, dist_params, test_fn=st.wald.ppf, test="icdf") From c124669a4621da403ddd113cef7c44f12f621be6 Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Mon, 20 Feb 2023 12:56:32 +0100 Subject: [PATCH 2/9] modified Wald.icdf to use only pytensor expressions, added docs --- pymc/distributions/continuous.py | 37 ++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 0e3baecc0f..86e6f40891 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -1019,8 +1019,41 @@ def logcdf(value, mu, lam, alpha): msg="mu > 0, lam > 0, alpha >= 0", ) - def icdf(mu, lam): - return stats.invgauss(mu=mu, scale=lam) + def icdf(value, mu, lam, alpha=0): + """Calculate the inverse cumulative distribution function (icdf) of the Wald distribution. + + Args + ---- + value (float): Probability value between 0 and 1. + mu (float): Mean of the distribution. + lam (float): Scale of the distribution. + + Returns + ------- + float: The value x such that P(W <= x) = p, where W is the Wald distribution with given mu and lam. + + Raises + ------ + ValueError: If value is not between 0 and 1. + ValueError: If lam is not positive. + """ + # Compute standard deviation and location parameter + std = at.sqrt(lam) + loc = alpha + mu + + # Compute inverse standard normal CDF + z = at.sqrt(2) * at.erfinv(2 * value - 1) + + # Compute Wald ICDF + x = loc + std * z / (1 - 0.5 * std * z) + + return check_parameters( + x, + 0 <= value <= 1, + lam > 0, + alpha >= 0, + msg="0<=val<=1, lam > 0, alpha >= 0", + ) class BetaClippedRV(BetaRV): From fa510fd6ab2174679425206254b43a4307c1f102 Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Mon, 20 Feb 2023 15:30:42 +0100 Subject: [PATCH 3/9] included suggestions from review --- pymc/distributions/continuous.py | 21 ++------------------- 1 file changed, 2 insertions(+), 19 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 86e6f40891..5595726590 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -1019,24 +1019,7 @@ def logcdf(value, mu, lam, alpha): msg="mu > 0, lam > 0, alpha >= 0", ) - def icdf(value, mu, lam, alpha=0): - """Calculate the inverse cumulative distribution function (icdf) of the Wald distribution. - - Args - ---- - value (float): Probability value between 0 and 1. - mu (float): Mean of the distribution. - lam (float): Scale of the distribution. - - Returns - ------- - float: The value x such that P(W <= x) = p, where W is the Wald distribution with given mu and lam. - - Raises - ------ - ValueError: If value is not between 0 and 1. - ValueError: If lam is not positive. - """ + def icdf(value, mu, lam, alpha): # Compute standard deviation and location parameter std = at.sqrt(lam) loc = alpha + mu @@ -1052,7 +1035,7 @@ def icdf(value, mu, lam, alpha=0): 0 <= value <= 1, lam > 0, alpha >= 0, - msg="0<=val<=1, lam > 0, alpha >= 0", + msg="0 <= val <= 1, lam > 0, alpha >= 0", ) From c113a09f06ac7def1466a2226deb2b31583335ad Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Mon, 20 Feb 2023 16:13:04 +0100 Subject: [PATCH 4/9] adapted previously created functions to match review suggestions, added Truncated Normal --- pymc/distributions/continuous.py | 25 +++++++++++++++++++-- pymc/tests/distributions/test_continuous.py | 8 +++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 5595726590..9b0af2f447 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -729,6 +729,27 @@ def logp(value, mu, sigma, lower, upper): return logp + def icdf(value, mu, sigma, lower, upper): + value = np.asarray(value) + eps = np.finfo(float).eps + + # Standardize the lower and upper bounds + lower_std = (lower - mu) / sigma + upper_std = (upper - mu) / sigma + + # Compute the cdf of the standard normal distribution at the bounds + Phi_a = 0.5 * (1 + np.math.erf(lower_std / np.sqrt(2))) + Phi_b = 0.5 * (1 + np.math.erf(upper_std / np.sqrt(2))) + + # Compute the cdf of the truncated normal distribution at the quantiles + Phi_q = Phi_a + value * (Phi_b - Phi_a) + + # Invert the cdf of the standard normal distribution at the quantiles + z = at.erfinv(2 * Phi_q - 1) * np.sqrt(2) + + # Transform the samples back to the truncated normal distribution + return mu + sigma * at.clip(z, lower_std + eps, upper_std - eps) + @_default_transform.register(TruncatedNormal) def truncated_normal_default_transform(op, rv): @@ -838,8 +859,8 @@ def logcdf(value, loc, sigma): msg="sigma > 0", ) - def icdf(value, lower, upper, loc, sigma): - return stats.truncnorm.ppf(q=value, a=lower, b=upper, loc=loc, scale=sigma) + def icdf(value, sigma): + return at.sqrt(2 * sigma**2) * at.erfinv(2 * value - 1) class WaldRV(RandomVariable): diff --git a/pymc/tests/distributions/test_continuous.py b/pymc/tests/distributions/test_continuous.py index 6745a65e75..8ffc5951c5 100644 --- a/pymc/tests/distributions/test_continuous.py +++ b/pymc/tests/distributions/test_continuous.py @@ -38,6 +38,7 @@ from pymc.logprob.abstract import logcdf from pymc.logprob.joint_logprob import logp from pymc.logprob.utils import ParameterValueError +from pymc.pymc.distributions.continuous import HalfNormal from pymc.pytensorf import floatX from pymc.tests.distributions.util import ( BaseTestDistributionRandom, @@ -2310,3 +2311,10 @@ def test_wald_icdf(self, *dist_params, obs, size): x = Wald.dist(*dist_params_at, size=size_at) scipy_logprob_tester(x, obs, dist_params, test_fn=st.wald.ppf, test="icdf") + + def test_half_normal_icdf(self, *dist_params, obs, size): + dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) + dist_params = dict(zip(dist_params_at, dist_params)) + + x = HalfNormal.dist(*dist_params_at, size=size_at) + scipy_logprob_tester(x, obs, dist_params, test_fn=st.wald.ppf, test="icdf") From 34673de586b4ae029d93239268baf6daf0b80a7b Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Tue, 21 Feb 2023 13:01:39 +0100 Subject: [PATCH 5/9] added test for uniform dist, work in progress on the rest --- pymc/distributions/continuous.py | 8 ++--- pymc/tests/distributions/test_continuous.py | 40 ++++++--------------- 2 files changed, 14 insertions(+), 34 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 9b0af2f447..058fa70b95 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -344,7 +344,7 @@ def logcdf(value, lower, upper): ) def icdf(value, lower, upper): - return lower + (upper - lower) * value + return (upper * value) + lower @_default_transform.register(Uniform) @@ -653,9 +653,9 @@ def dist( cls, mu: Optional[DIST_PARAMETER_TYPES] = None, sigma: Optional[DIST_PARAMETER_TYPES] = None, - tau: Optional[DIST_PARAMETER_TYPES] = None, lower: Optional[DIST_PARAMETER_TYPES] = None, upper: Optional[DIST_PARAMETER_TYPES] = None, + tau: Optional[DIST_PARAMETER_TYPES] = None, *args, **kwargs, ) -> RandomVariable: @@ -738,8 +738,8 @@ def icdf(value, mu, sigma, lower, upper): upper_std = (upper - mu) / sigma # Compute the cdf of the standard normal distribution at the bounds - Phi_a = 0.5 * (1 + np.math.erf(lower_std / np.sqrt(2))) - Phi_b = 0.5 * (1 + np.math.erf(upper_std / np.sqrt(2))) + Phi_a = 0.5 * (1 + at.erf(lower_std / at.sqrt(2))) + Phi_b = 0.5 * (1 + at.erf(upper_std / at.sqrt(2))) # Compute the cdf of the truncated normal distribution at the quantiles Phi_q = Phi_a + value * (Phi_b - Phi_a) diff --git a/pymc/tests/distributions/test_continuous.py b/pymc/tests/distributions/test_continuous.py index 8ffc5951c5..9863d9b0bc 100644 --- a/pymc/tests/distributions/test_continuous.py +++ b/pymc/tests/distributions/test_continuous.py @@ -26,19 +26,11 @@ import pymc as pm -from pymc.distributions.continuous import ( - Normal, - TruncatedNormal, - Uniform, - Wald, - get_tau_sigma, - interpolated, -) +from pymc.distributions.continuous import Normal, Uniform, get_tau_sigma, interpolated from pymc.distributions.dist_math import clipped_beta_rvs from pymc.logprob.abstract import logcdf from pymc.logprob.joint_logprob import logp from pymc.logprob.utils import ParameterValueError -from pymc.pymc.distributions.continuous import HalfNormal from pymc.pytensorf import floatX from pymc.tests.distributions.util import ( BaseTestDistributionRandom, @@ -2289,32 +2281,20 @@ def test_normal_icdf(self, dist_params, obs, size): dist_params = dict(zip(dist_params_at, dist_params)) x = Normal.dist(*dist_params_at, size=size_at) + scipy_logprob_tester(x, obs, dist_params, test_fn=st.norm.ppf, test="icdf") + @pytest.mark.parametrize( + "dist_params, obs, size", + [ + ((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), + ((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), (4, 6)), + ((0, 10), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), + ], + ) def test_uniform_icdf(self, dist_params, obs, size): dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) dist_params = dict(zip(dist_params_at, dist_params)) x = Uniform.dist(*dist_params_at, size=size_at) scipy_logprob_tester(x, obs, dist_params, test_fn=st.uniform.ppf, test="icdf") - - def test_truncated_normal_icdf(self, dist_params, obs, size): - dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) - dist_params = dict(zip(dist_params_at, dist_params)) - - x = TruncatedNormal.dist(*dist_params_at, size=size_at) - scipy_logprob_tester(x, obs, dist_params, test_fn=st.truncnorm.ppf, test="icdf") - - def test_wald_icdf(self, *dist_params, obs, size): - dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) - dist_params = dict(zip(dist_params_at, dist_params)) - - x = Wald.dist(*dist_params_at, size=size_at) - scipy_logprob_tester(x, obs, dist_params, test_fn=st.wald.ppf, test="icdf") - - def test_half_normal_icdf(self, *dist_params, obs, size): - dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) - dist_params = dict(zip(dist_params_at, dist_params)) - - x = HalfNormal.dist(*dist_params_at, size=size_at) - scipy_logprob_tester(x, obs, dist_params, test_fn=st.wald.ppf, test="icdf") From d87dc9550b911e02a3a7c249dd4b3f0c2a78fb78 Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Mon, 6 Mar 2023 10:36:24 +0100 Subject: [PATCH 6/9] added tests for uniform icdf --- pymc/distributions/continuous.py | 8 +++---- pymc/tests/distributions/test_continuous.py | 26 ++++++++++++++++++--- 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 058fa70b95..102484c350 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -344,7 +344,7 @@ def logcdf(value, lower, upper): ) def icdf(value, lower, upper): - return (upper * value) + lower + return lower + (upper - lower) * value @_default_transform.register(Uniform) @@ -653,9 +653,9 @@ def dist( cls, mu: Optional[DIST_PARAMETER_TYPES] = None, sigma: Optional[DIST_PARAMETER_TYPES] = None, + tau: Optional[DIST_PARAMETER_TYPES] = None, lower: Optional[DIST_PARAMETER_TYPES] = None, upper: Optional[DIST_PARAMETER_TYPES] = None, - tau: Optional[DIST_PARAMETER_TYPES] = None, *args, **kwargs, ) -> RandomVariable: @@ -859,8 +859,8 @@ def logcdf(value, loc, sigma): msg="sigma > 0", ) - def icdf(value, sigma): - return at.sqrt(2 * sigma**2) * at.erfinv(2 * value - 1) + def icdf(value, loc, sigma): + return loc + sigma * np.sqrt(2) * np.erfinv(2 * value - 1) class WaldRV(RandomVariable): diff --git a/pymc/tests/distributions/test_continuous.py b/pymc/tests/distributions/test_continuous.py index 9863d9b0bc..aa17dc502c 100644 --- a/pymc/tests/distributions/test_continuous.py +++ b/pymc/tests/distributions/test_continuous.py @@ -2288,13 +2288,33 @@ def test_normal_icdf(self, dist_params, obs, size): "dist_params, obs, size", [ ((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), - ((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), (4, 6)), + ((-1, 2), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), ((0, 10), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), ], ) def test_uniform_icdf(self, dist_params, obs, size): dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) dist_params = dict(zip(dist_params_at, dist_params)) + x = Uniform.dist(*dist_params_at) + scipy_logprob_tester( + x, + obs, + dist_params, + test_fn=lambda val, lower, upper: st.uniform.ppf(val, lower, upper - lower), + test="icdf", + ) + + """@pytest.mark.parametrize( + "dist_params, obs, size", + [ + ((0, 1), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), + #((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), (4, 6)), + #((0, 10), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), + ], + ) + def test_half_normal_icdf(self, *dist_params, obs, size): + dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) + dist_params = dict(zip(dist_params_at, dist_params)) - x = Uniform.dist(*dist_params_at, size=size_at) - scipy_logprob_tester(x, obs, dist_params, test_fn=st.uniform.ppf, test="icdf") + x = HalfNormal.dist(*dist_params_at, size=size_at) + scipy_logprob_tester(x, obs, dist_params, test_fn=st.truncnorm.ppf, test="icdf")""" From 4d636560d7d3d8dee3a4b4f5420f26b2735f0640 Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Thu, 9 Mar 2023 11:31:46 +0100 Subject: [PATCH 7/9] removing icdf's that don't pass tests yet --- pymc/distributions/continuous.py | 22 --------------------- pymc/tests/distributions/test_continuous.py | 15 -------------- 2 files changed, 37 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 102484c350..266ab51a00 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -859,9 +859,6 @@ def logcdf(value, loc, sigma): msg="sigma > 0", ) - def icdf(value, loc, sigma): - return loc + sigma * np.sqrt(2) * np.erfinv(2 * value - 1) - class WaldRV(RandomVariable): name = "wald" @@ -1040,25 +1037,6 @@ def logcdf(value, mu, lam, alpha): msg="mu > 0, lam > 0, alpha >= 0", ) - def icdf(value, mu, lam, alpha): - # Compute standard deviation and location parameter - std = at.sqrt(lam) - loc = alpha + mu - - # Compute inverse standard normal CDF - z = at.sqrt(2) * at.erfinv(2 * value - 1) - - # Compute Wald ICDF - x = loc + std * z / (1 - 0.5 * std * z) - - return check_parameters( - x, - 0 <= value <= 1, - lam > 0, - alpha >= 0, - msg="0 <= val <= 1, lam > 0, alpha >= 0", - ) - class BetaClippedRV(BetaRV): @classmethod diff --git a/pymc/tests/distributions/test_continuous.py b/pymc/tests/distributions/test_continuous.py index aa17dc502c..6e0ed5e28b 100644 --- a/pymc/tests/distributions/test_continuous.py +++ b/pymc/tests/distributions/test_continuous.py @@ -2303,18 +2303,3 @@ def test_uniform_icdf(self, dist_params, obs, size): test_fn=lambda val, lower, upper: st.uniform.ppf(val, lower, upper - lower), test="icdf", ) - - """@pytest.mark.parametrize( - "dist_params, obs, size", - [ - ((0, 1), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), - #((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), (4, 6)), - #((0, 10), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), - ], - ) - def test_half_normal_icdf(self, *dist_params, obs, size): - dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) - dist_params = dict(zip(dist_params_at, dist_params)) - - x = HalfNormal.dist(*dist_params_at, size=size_at) - scipy_logprob_tester(x, obs, dist_params, test_fn=st.truncnorm.ppf, test="icdf")""" From ecbb89a47e95b8e36bb3710cb8a43c7de1ae0d0b Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Sat, 18 Mar 2023 15:18:32 +0300 Subject: [PATCH 8/9] adjusting uniform icdf to fit the changes introduced in #6583 --- pymc/distributions/continuous.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 2d27f02347..f970aacb10 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -346,7 +346,9 @@ def logcdf(value, lower, upper): ) def icdf(value, lower, upper): - return lower + (upper - lower) * value + res = lower + (upper - lower) * value + res = check_icdf_value(res, value) + return check_icdf_parameters(res, lower < upper) @_default_transform.register(Uniform) From 63641a79ac04aa94bffd788c7a355bb2b5aa6df7 Mon Sep 17 00:00:00 2001 From: Michal Raczycki Date: Sun, 19 Mar 2023 12:25:30 +0300 Subject: [PATCH 9/9] adapted uniform.icdf tests to run on new test suite, extended testing.py check_icdf with skip_paradomain_outside_edge_test param --- pymc/distributions/continuous.py | 21 ----------- pymc/testing.py | 31 ++++++++++------- tests/distributions/test_continuous.py | 48 +++++--------------------- 3 files changed, 27 insertions(+), 73 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index f970aacb10..b45eaeb96f 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -739,27 +739,6 @@ def logp(value, mu, sigma, lower, upper): return logp - def icdf(value, mu, sigma, lower, upper): - value = np.asarray(value) - eps = np.finfo(float).eps - - # Standardize the lower and upper bounds - lower_std = (lower - mu) / sigma - upper_std = (upper - mu) / sigma - - # Compute the cdf of the standard normal distribution at the bounds - Phi_a = 0.5 * (1 + at.erf(lower_std / at.sqrt(2))) - Phi_b = 0.5 * (1 + at.erf(upper_std / at.sqrt(2))) - - # Compute the cdf of the truncated normal distribution at the quantiles - Phi_q = Phi_a + value * (Phi_b - Phi_a) - - # Invert the cdf of the standard normal distribution at the quantiles - z = at.erfinv(2 * Phi_q - 1) * np.sqrt(2) - - # Transform the samples back to the truncated normal distribution - return mu + sigma * at.clip(z, lower_std + eps, upper_std - eps) - @_default_transform.register(TruncatedNormal) def truncated_normal_default_transform(op, rv): diff --git a/pymc/testing.py b/pymc/testing.py index ea3ccfc46f..3bb222222f 100644 --- a/pymc/testing.py +++ b/pymc/testing.py @@ -526,6 +526,7 @@ def check_icdf( pymc_dist: Distribution, paramdomains: Dict[str, Domain], scipy_icdf: Callable, + skip_paramdomain_outside_edge_test=False, decimal: Optional[int] = None, n_samples: int = 100, ) -> None: @@ -548,7 +549,7 @@ def check_icdf( paramdomains : Dictionary of Parameter : Domain pairs Supported domains of distribution parameters scipy_icdf : Scipy icdf method - Scipy icdf (ppp) method of equivalent pymc_dist distribution + Scipy icdf (ppf) method of equivalent pymc_dist distribution decimal : int, optional Level of precision with which pymc_dist and scipy_icdf are compared. Defaults to 6 for float64 and 3 for float32 @@ -557,6 +558,9 @@ def check_icdf( are compared between pymc and scipy methods. If n_samples is below the total number of combinations, a random subset is evaluated. Setting n_samples = -1, will return all possible combinations. Defaults to 100 + skip_paradomain_outside_edge_test : Bool + Whether to run test 2., which checks that pymc distribution icdf + returns nan for invalid parameter values outside the supported domain edge """ if decimal is None: @@ -586,19 +590,20 @@ def check_icdf( valid_params = {param: paramdomain.vals[0] for param, paramdomain in paramdomains.items()} valid_params["q"] = valid_value - # Test pymc distribution raises ParameterValueError for parameters outside the - # supported domain edges (excluding edges) - invalid_params = find_invalid_scalar_params(paramdomains) - for invalid_param, invalid_edges in invalid_params.items(): - for invalid_edge in invalid_edges: - if invalid_edge is None: - continue + if not skip_paramdomain_outside_edge_test: + # Test pymc distribution raises ParameterValueError for parameters outside the + # supported domain edges (excluding edges) + invalid_params = find_invalid_scalar_params(paramdomains) + for invalid_param, invalid_edges in invalid_params.items(): + for invalid_edge in invalid_edges: + if invalid_edge is None: + continue - point = valid_params.copy() - point[invalid_param] = invalid_edge - with pytest.raises(ParameterValueError): - pymc_icdf(**point) - pytest.fail(f"test_params={point}") + point = valid_params.copy() + point[invalid_param] = invalid_edge + with pytest.raises(ParameterValueError): + pymc_icdf(**point) + pytest.fail(f"test_params={point}") # Test that values below 0 or above 1 evaluate to nan invalid_values = find_invalid_scalar_params({"q": domain})["q"] diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index ea84afef91..8b4484a66c 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -28,7 +28,7 @@ from pymc.distributions.continuous import Normal, Uniform, get_tau_sigma, interpolated from pymc.distributions.dist_math import clipped_beta_rvs -from pymc.logprob.abstract import logcdf +from pymc.logprob.abstract import icdf, logcdf from pymc.logprob.joint_logprob import logp from pymc.logprob.utils import ParameterValueError from pymc.pytensorf import floatX @@ -176,6 +176,12 @@ def test_uniform(self): lambda value, lower, upper: st.uniform.logcdf(value, lower, upper - lower), skip_paramdomain_outside_edge_test=True, ) + check_icdf( + pm.Uniform, + {"lower": -Rplusunif, "upper": Rplusunif}, + lambda q, lower, upper: st.uniform.ppf(q=q, loc=lower, scale=upper - lower), + skip_paramdomain_outside_edge_test=True, + ) # Custom logp / logcdf check for invalid parameters invalid_dist = pm.Uniform.dist(lower=1, upper=0) with pytensor.config.change_flags(mode=Mode("py")): @@ -183,6 +189,8 @@ def test_uniform(self): logp(invalid_dist, np.array(0.5)).eval() with pytest.raises(ParameterValueError): logcdf(invalid_dist, np.array(0.5)).eval() + with pytest.raises(ParameterValueError): + icdf(invalid_dist, np.array(0.5)).eval() def test_triangular(self): check_logp( @@ -2277,41 +2285,3 @@ def dist(cls, **kwargs): extra_args={"rng": pytensor.shared(rng)}, ref_rand=ref_rand, ) - - -class TestICDF: - @pytest.mark.parametrize( - "dist_params, obs, size", - [ - ((0, 1), np.array([-0.5, 0, 0.3, 0.5, 1, 1.5], dtype=np.float64), ()), - ((-1, 20), np.array([-0.5, 0, 0.3, 0.5, 1, 1.5], dtype=np.float64), ()), - ((-1, 20), np.array([-0.5, 0, 0.3, 0.5, 1, 1.5], dtype=np.float64), (2, 3)), - ], - ) - def test_normal_icdf(self, dist_params, obs, size): - dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) - dist_params = dict(zip(dist_params_at, dist_params)) - - x = Normal.dist(*dist_params_at, size=size_at) - - scipy_logprob_tester(x, obs, dist_params, test_fn=st.norm.ppf, test="icdf") - - @pytest.mark.parametrize( - "dist_params, obs, size", - [ - ((-5, 4), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), - ((-1, 2), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), - ((0, 10), np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64), ()), - ], - ) - def test_uniform_icdf(self, dist_params, obs, size): - dist_params_at, obs_at, size_at = create_pytensor_params(dist_params, obs, size) - dist_params = dict(zip(dist_params_at, dist_params)) - x = Uniform.dist(*dist_params_at) - scipy_logprob_tester( - x, - obs, - dist_params, - test_fn=lambda val, lower, upper: st.uniform.ppf(val, lower, upper - lower), - test="icdf", - )