From 9ecb0083b964430e83b2754d49db450a3f9b3592 Mon Sep 17 00:00:00 2001 From: Gokul D Date: Tue, 4 Apr 2023 01:09:28 +0530 Subject: [PATCH 1/6] Added ICDF for the Kumaraswamy distribution. --- pymc/distributions/continuous.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index fd6e414605..c0c957f20e 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -1323,6 +1323,16 @@ def logcdf(value, a, b): msg="a > 0, b > 0", ) + def icdf(value, a, b): + res = pt.exp(pt.reciprocal(a) * pt.log1mexp(pt.reciprocal(b) * pt.log1p(-value))) + res = check_icdf_value(res, value) + return check_icdf_parameters( + res, + a > 0, + b > 0, + msg="a > 0, b > 0", + ) + class Exponential(PositiveContinuous): r""" From edd7cb66f19bf1c7bd9e831d590946df7ba59571 Mon Sep 17 00:00:00 2001 From: Gokul D Date: Tue, 4 Apr 2023 01:13:18 +0530 Subject: [PATCH 2/6] Added icdf - logcdf consistency tests. --- pymc/testing.py | 76 ++++++++++++++++++++++++++ tests/distributions/test_continuous.py | 5 ++ tests/distributions/test_discrete.py | 11 ++++ 3 files changed, 92 insertions(+) diff --git a/pymc/testing.py b/pymc/testing.py index 0405cc2f6c..8df1494182 100644 --- a/pymc/testing.py +++ b/pymc/testing.py @@ -668,6 +668,82 @@ def check_selfconsistency_discrete_logcdf( ) +def check_selfconsistency_continuous_icdf( + distribution: Distribution, + paramdomains: Dict[str, Domain], + decimal: Optional[int] = None, + n_samples: int = 100, +) -> None: + """ + Check that the icdf and logcdf functions of the distribution are consistent for a sample of probability values. + """ + if decimal is None: + decimal = select_by_precision(float64=6, float32=3) + + dist = create_dist_from_paramdomains(distribution, paramdomains) + value = dist.type() + value.name = "value" + + dist_icdf = icdf(dist, value) + dist_icdf_fn = pytensor.function(list(inputvars(dist_icdf)), dist_icdf) + + dist_logcdf = logcdf(dist, value) + dist_logcdf_fn = compile_pymc(list(inputvars(dist_logcdf)), dist_logcdf) + + domains = paramdomains.copy() + domains["value"] = Domain(np.linspace(0, 1, 10)) + + for point in product(domains, n_samples=n_samples): + point = dict(point) + value = point.pop("value") + + with pytensor.config.change_flags(mode=Mode("py")): + npt.assert_almost_equal( + value, + np.exp(dist_logcdf_fn(**point, value=dist_icdf_fn(**point, value=value))), + decimal=decimal, + err_msg=f"point: {point}, value: {value}", + ) + + +def check_selfconsistency_discrete_icdf( + distribution: Distribution, + domain: Domain, + paramdomains: Dict[str, Domain], + n_samples: int = 100, +) -> None: + """ + Check that the icdf and logcdf functions of the distribution are + consistent for a sample of values in the domain of the + distribution. + """ + dist = create_dist_from_paramdomains(distribution, paramdomains) + + value = pt.TensorType(dtype="float64", shape=[])("value") + + dist_icdf = icdf(dist, value) + dist_icdf_fn = pytensor.function(list(inputvars(dist_icdf)), dist_icdf) + + dist_logcdf = logcdf(dist, value) + dist_logcdf_fn = compile_pymc(list(inputvars(dist_logcdf)), dist_logcdf) + + domains = paramdomains.copy() + domains["value"] = domain + + for point in product(domains, n_samples=n_samples): + point = dict(point) + value = point.pop("value") + + with pytensor.config.change_flags(mode=Mode("py")): + expected_value = value + computed_value = dist_icdf_fn( + **point, value=np.exp(dist_logcdf_fn(**point, value=value)) + ) + assert ( + expected_value == computed_value + ), f"expected_value = {expected_value}, computed_value = {computed_value}, {point}" + + def assert_support_point_is_expected(model, expected, check_finite_logp=True): fn = make_initial_point_fn( model=model, diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index 7209382666..35def0ff6a 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -45,6 +45,7 @@ check_icdf, check_logcdf, check_logp, + check_selfconsistency_continuous_icdf, continuous_random_tester, seeded_numpy_distribution_builder, seeded_scipy_distribution_builder, @@ -441,6 +442,10 @@ def scipy_log_cdf(value, a, b): {"a": Rplus, "b": Rplus}, scipy_log_cdf, ) + check_selfconsistency_continuous_icdf( + pm.Kumaraswamy, + {"a": Rplusbig, "b": Rplusbig}, + ) def test_exponential(self): check_logp( diff --git a/tests/distributions/test_discrete.py b/tests/distributions/test_discrete.py index 55e8c23128..06c3d2e75e 100644 --- a/tests/distributions/test_discrete.py +++ b/tests/distributions/test_discrete.py @@ -53,6 +53,7 @@ check_icdf, check_logcdf, check_logp, + check_selfconsistency_discrete_icdf, check_selfconsistency_discrete_logcdf, seeded_numpy_distribution_builder, seeded_scipy_distribution_builder, @@ -121,6 +122,11 @@ def test_discrete_unif(self): lambda q, lower, upper: st.randint.ppf(q=q, low=lower, high=upper + 1), skip_paramdomain_outside_edge_test=True, ) + check_selfconsistency_discrete_icdf( + pm.DiscreteUniform, + Rdunif, + {"lower": -Rplusdunif, "upper": Rplusdunif}, + ) # Custom logp / logcdf check for invalid parameters invalid_dist = pm.DiscreteUniform.dist(lower=1, upper=0) with pytensor.config.change_flags(mode=Mode("py")): @@ -154,6 +160,11 @@ def test_geometric(self): {"p": Unit}, st.geom.ppf, ) + check_selfconsistency_discrete_icdf( + pm.Geometric, + Nat, + {"p": Unit}, + ) def test_hypergeometric(self): def modified_scipy_hypergeom_logcdf(value, N, k, n): From 603afad847bebd8f853066aa772c9d92b4874cb3 Mon Sep 17 00:00:00 2001 From: Gokul D Date: Sat, 29 Apr 2023 17:51:23 +0530 Subject: [PATCH 3/6] Tmp testing with np assert_almost_equal print messages. --- pymc/testing.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pymc/testing.py b/pymc/testing.py index 8df1494182..59d584e8fd 100644 --- a/pymc/testing.py +++ b/pymc/testing.py @@ -739,9 +739,11 @@ def check_selfconsistency_discrete_icdf( computed_value = dist_icdf_fn( **point, value=np.exp(dist_logcdf_fn(**point, value=value)) ) - assert ( - expected_value == computed_value - ), f"expected_value = {expected_value}, computed_value = {computed_value}, {point}" + npt.assert_almost_equal( + expected_value, + computed_value, + err_msg=f"point: {point}, value: {value}, cdf: {np.exp(dist_logcdf_fn(**point, value=value))}", + ) def assert_support_point_is_expected(model, expected, check_finite_logp=True): From 64b9bc7405c45f1b8cbcedf48451b9086de3b013 Mon Sep 17 00:00:00 2001 From: Gokul D Date: Sat, 29 Apr 2023 22:23:26 +0530 Subject: [PATCH 4/6] Testing a fix by truncating cdf to precision. --- pymc/testing.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pymc/testing.py b/pymc/testing.py index 59d584e8fd..04db416fe2 100644 --- a/pymc/testing.py +++ b/pymc/testing.py @@ -710,6 +710,7 @@ def check_selfconsistency_discrete_icdf( distribution: Distribution, domain: Domain, paramdomains: Dict[str, Domain], + decimal: Optional[int] = None, n_samples: int = 100, ) -> None: """ @@ -717,6 +718,13 @@ def check_selfconsistency_discrete_icdf( consistent for a sample of values in the domain of the distribution. """ + + def ftrunc(values, decimal=0): + return np.trunc(values * 10**decimal) / (10**decimal) + + if decimal is None: + decimal = select_by_precision(float64=6, float32=3) + dist = create_dist_from_paramdomains(distribution, paramdomains) value = pt.TensorType(dtype="float64", shape=[])("value") @@ -737,7 +745,7 @@ def check_selfconsistency_discrete_icdf( with pytensor.config.change_flags(mode=Mode("py")): expected_value = value computed_value = dist_icdf_fn( - **point, value=np.exp(dist_logcdf_fn(**point, value=value)) + **point, value=ftrunc(np.exp(dist_logcdf_fn(**point, value=value)), decimal=decimal) ) npt.assert_almost_equal( expected_value, From 6a7111a97f4fea93266a6fc8c464cf6053715eb0 Mon Sep 17 00:00:00 2001 From: Gokul D Date: Thu, 1 Jun 2023 14:23:09 +0530 Subject: [PATCH 5/6] Revert "Tmp testing with np assert_almost_equal print messages." This reverts commit b0c93b77e3d39c81b22325b644cf94c2beeb9e9c. Avoid mypy test failure. This is not needed, we have identified the root cause now. --- pymc/testing.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/pymc/testing.py b/pymc/testing.py index 04db416fe2..cb64170de3 100644 --- a/pymc/testing.py +++ b/pymc/testing.py @@ -747,11 +747,9 @@ def ftrunc(values, decimal=0): computed_value = dist_icdf_fn( **point, value=ftrunc(np.exp(dist_logcdf_fn(**point, value=value)), decimal=decimal) ) - npt.assert_almost_equal( - expected_value, - computed_value, - err_msg=f"point: {point}, value: {value}, cdf: {np.exp(dist_logcdf_fn(**point, value=value))}", - ) + assert ( + expected_value == computed_value + ), f"expected_value = {expected_value}, computed_value = {computed_value}, {point}" def assert_support_point_is_expected(model, expected, check_finite_logp=True): From 813aece21a10fc6980b8d1b7f512df78f85c3ca1 Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Sat, 8 Nov 2025 13:48:19 +0100 Subject: [PATCH 6/6] Remove discrete selfconsistency icdf --- pymc/testing.py | 92 ++++++++------------------ tests/distributions/test_continuous.py | 4 +- tests/distributions/test_discrete.py | 11 --- 3 files changed, 29 insertions(+), 78 deletions(-) diff --git a/pymc/testing.py b/pymc/testing.py index cb64170de3..551356ee33 100644 --- a/pymc/testing.py +++ b/pymc/testing.py @@ -668,88 +668,50 @@ def check_selfconsistency_discrete_logcdf( ) -def check_selfconsistency_continuous_icdf( +def check_selfconsistency_icdf( distribution: Distribution, - paramdomains: Dict[str, Domain], - decimal: Optional[int] = None, + paramdomains: dict[str, Domain], + *, + decimal: int | None = None, n_samples: int = 100, ) -> None: - """ - Check that the icdf and logcdf functions of the distribution are consistent for a sample of probability values. - """ - if decimal is None: - decimal = select_by_precision(float64=6, float32=3) - - dist = create_dist_from_paramdomains(distribution, paramdomains) - value = dist.type() - value.name = "value" - - dist_icdf = icdf(dist, value) - dist_icdf_fn = pytensor.function(list(inputvars(dist_icdf)), dist_icdf) - - dist_logcdf = logcdf(dist, value) - dist_logcdf_fn = compile_pymc(list(inputvars(dist_logcdf)), dist_logcdf) - - domains = paramdomains.copy() - domains["value"] = Domain(np.linspace(0, 1, 10)) - - for point in product(domains, n_samples=n_samples): - point = dict(point) - value = point.pop("value") - - with pytensor.config.change_flags(mode=Mode("py")): - npt.assert_almost_equal( - value, - np.exp(dist_logcdf_fn(**point, value=dist_icdf_fn(**point, value=value))), - decimal=decimal, - err_msg=f"point: {point}, value: {value}", - ) + """Check that the icdf and logcdf functions of the distribution are consistent. - -def check_selfconsistency_discrete_icdf( - distribution: Distribution, - domain: Domain, - paramdomains: Dict[str, Domain], - decimal: Optional[int] = None, - n_samples: int = 100, -) -> None: - """ - Check that the icdf and logcdf functions of the distribution are - consistent for a sample of values in the domain of the - distribution. + Only works with continuous distributions. """ - - def ftrunc(values, decimal=0): - return np.trunc(values * 10**decimal) / (10**decimal) - if decimal is None: decimal = select_by_precision(float64=6, float32=3) dist = create_dist_from_paramdomains(distribution, paramdomains) - - value = pt.TensorType(dtype="float64", shape=[])("value") - + if dist.type.dtype.startswith("int"): + raise NotImplementedError( + "check_selfconsistency_icdf is not robust against discrete distributions." + ) + value = dist.astype("float64").type("value") dist_icdf = icdf(dist, value) - dist_icdf_fn = pytensor.function(list(inputvars(dist_icdf)), dist_icdf) + dist_cdf = pt.exp(logcdf(dist, value)) - dist_logcdf = logcdf(dist, value) - dist_logcdf_fn = compile_pymc(list(inputvars(dist_logcdf)), dist_logcdf) + py_mode = Mode("py") + dist_icdf_fn = pytensor.function(list(inputvars(dist_icdf)), dist_icdf, mode=py_mode) + dist_cdf_fn = compile(list(inputvars(dist_cdf)), dist_cdf, mode=py_mode) domains = paramdomains.copy() - domains["value"] = domain + domains["value"] = Domain(np.linspace(0, 1, 10)) for point in product(domains, n_samples=n_samples): point = dict(point) value = point.pop("value") - - with pytensor.config.change_flags(mode=Mode("py")): - expected_value = value - computed_value = dist_icdf_fn( - **point, value=ftrunc(np.exp(dist_logcdf_fn(**point, value=value)), decimal=decimal) - ) - assert ( - expected_value == computed_value - ), f"expected_value = {expected_value}, computed_value = {computed_value}, {point}" + icdf_value = dist_icdf_fn(**point, value=value) + recovered_value = dist_cdf_fn( + **point, + value=icdf_value, + ) + np.testing.assert_almost_equal( + value, + recovered_value, + decimal=decimal, + err_msg=f"point: {point}", + ) def assert_support_point_is_expected(model, expected, check_finite_logp=True): diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index 35def0ff6a..ef21aea1a3 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -45,7 +45,7 @@ check_icdf, check_logcdf, check_logp, - check_selfconsistency_continuous_icdf, + check_selfconsistency_icdf, continuous_random_tester, seeded_numpy_distribution_builder, seeded_scipy_distribution_builder, @@ -442,7 +442,7 @@ def scipy_log_cdf(value, a, b): {"a": Rplus, "b": Rplus}, scipy_log_cdf, ) - check_selfconsistency_continuous_icdf( + check_selfconsistency_icdf( pm.Kumaraswamy, {"a": Rplusbig, "b": Rplusbig}, ) diff --git a/tests/distributions/test_discrete.py b/tests/distributions/test_discrete.py index 06c3d2e75e..55e8c23128 100644 --- a/tests/distributions/test_discrete.py +++ b/tests/distributions/test_discrete.py @@ -53,7 +53,6 @@ check_icdf, check_logcdf, check_logp, - check_selfconsistency_discrete_icdf, check_selfconsistency_discrete_logcdf, seeded_numpy_distribution_builder, seeded_scipy_distribution_builder, @@ -122,11 +121,6 @@ def test_discrete_unif(self): lambda q, lower, upper: st.randint.ppf(q=q, low=lower, high=upper + 1), skip_paramdomain_outside_edge_test=True, ) - check_selfconsistency_discrete_icdf( - pm.DiscreteUniform, - Rdunif, - {"lower": -Rplusdunif, "upper": Rplusdunif}, - ) # Custom logp / logcdf check for invalid parameters invalid_dist = pm.DiscreteUniform.dist(lower=1, upper=0) with pytensor.config.change_flags(mode=Mode("py")): @@ -160,11 +154,6 @@ def test_geometric(self): {"p": Unit}, st.geom.ppf, ) - check_selfconsistency_discrete_icdf( - pm.Geometric, - Nat, - {"p": Unit}, - ) def test_hypergeometric(self): def modified_scipy_hypergeom_logcdf(value, N, k, n):