From 711ee4b7b04afa554353bd4727b4b75acddf945d Mon Sep 17 00:00:00 2001 From: amyoshino Date: Fri, 9 Jun 2023 15:51:44 -0400 Subject: [PATCH 1/7] adding icdf function and tests --- pymc/distributions/continuous.py | 9 +++++++++ tests/distributions/test_continuous.py | 10 ++++++++++ 2 files changed, 19 insertions(+) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index d597299997..56436edf02 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -1721,6 +1721,15 @@ def logcdf(value, mu, sigma): msg="sigma > 0", ) + def icdf(value, mu, sigma): + res = np.exp(mu + sigma * -np.sqrt(2.0) * pt.erfcinv(2 * value)) + res = check_icdf_value(res, value) + return check_icdf_parameters( + res, + sigma > 0, + msg="sigma > 0", + ) + Lognormal = LogNormal diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index b84a77e049..3b233cc85f 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -502,6 +502,16 @@ def test_lognormal(self): {"mu": R, "sigma": Rplusbig}, lambda value, mu, sigma: st.lognorm.logcdf(value, sigma, 0, np.exp(mu)), ) + check_icdf( + pm.LogNormal, + {"mu": R, "tau": Rplusbig}, + lambda q, mu, tau: floatX(st.lognorm.ppf(q, tau**-0.5, 0, np.exp(mu))), + ) + check_icdf( + pm.LogNormal, + {"mu": R, "sigma": Rplusbig}, + lambda q, mu, sigma: floatX(st.lognorm.ppf(q, sigma, 0, np.exp(mu))), + ) def test_studentt_logp(self): check_logp( From c9e48ec43d7d7968f79e4ccc9c071cd6eebb7de3 Mon Sep 17 00:00:00 2001 From: amyoshino Date: Mon, 19 Jun 2023 16:45:47 -0400 Subject: [PATCH 2/7] adding adjustments on tests for the lognormal icdf function --- pymc/distributions/continuous.py | 8 ++++---- tests/distributions/test_continuous.py | 7 ++++++- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 47d2a9d73d..ee70975b82 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -1714,6 +1714,7 @@ def logcdf(value, mu, sigma): -np.inf, normal_lcdf(mu, sigma, pt.log(value)), ) + return check_parameters( res, sigma > 0, @@ -1721,7 +1722,7 @@ def logcdf(value, mu, sigma): ) def icdf(value, mu, sigma): - res = np.exp(mu + sigma * -np.sqrt(2.0) * pt.erfcinv(2 * value)) + res = pt.exp(Normal.icdf(value, mu, sigma)) res = check_icdf_value(res, value) return check_icdf_parameters( res, @@ -3823,8 +3824,7 @@ class PolyaGamma(PositiveContinuous): import matplotlib.pyplot as plt import numpy as np from polyagamma import polyagamma_pdf - import arviz as az - plt.style.use('arviz-darkgrid') + plt.style.use('seaborn-darkgrid') x = np.linspace(0.01, 5, 500);x.sort() hs = [1., 5., 10., 15.] zs = [0.] * 4 @@ -3838,7 +3838,7 @@ class PolyaGamma(PositiveContinuous): ======== ============================= Support :math:`x \in (0, \infty)` - Mean :math:`\dfrac{h}{4}` if :math:`z=0`, :math:`\dfrac{tanh(z/2)h}{2z}` otherwise. + Mean :math:`dfrac{h}{4} if :math:`z=0`, :math:`\dfrac{tanh(z/2)h}{2z}` otherwise. Variance :math:`0.041666688h` if :math:`z=0`, :math:`\dfrac{h(sinh(z) - z)(1 - tanh^2(z/2))}{4z^3}` otherwise. ======== ============================= diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index 26c6ce08fb..3c5dd5fe4c 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -507,10 +507,15 @@ def test_lognormal(self): {"mu": R, "tau": Rplusbig}, lambda q, mu, tau: floatX(st.lognorm.ppf(q, tau**-0.5, 0, np.exp(mu))), ) + # Because we exponentiate the normal quantile function, setting sigma >= 9.5 + # return extreme values that results in relative errors above 4 digits + # we circumvent it by keeping it below or equal to 9. + custom_rplusbig = Domain([0, 0.5, 0.9, 0.99, 1, 1.5, 2, 9, np.inf]) check_icdf( pm.LogNormal, - {"mu": R, "sigma": Rplusbig}, + {"mu": R, "sigma": custom_rplusbig}, lambda q, mu, sigma: floatX(st.lognorm.ppf(q, sigma, 0, np.exp(mu))), + decimal=select_by_precision(float64=4, float32=3), ) def test_studentt_logp(self): From 1035cc6d389689824ae747e48e3cc00374fd1276 Mon Sep 17 00:00:00 2001 From: amyoshino Date: Mon, 19 Jun 2023 16:53:09 -0400 Subject: [PATCH 3/7] fixing changes wrong --- pymc/distributions/continuous.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index ee70975b82..11ce4ae180 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -3824,7 +3824,8 @@ class PolyaGamma(PositiveContinuous): import matplotlib.pyplot as plt import numpy as np from polyagamma import polyagamma_pdf - plt.style.use('seaborn-darkgrid') + import arviz as az + plt.style.use('arviz-darkgrid') x = np.linspace(0.01, 5, 500);x.sort() hs = [1., 5., 10., 15.] zs = [0.] * 4 @@ -3838,7 +3839,7 @@ class PolyaGamma(PositiveContinuous): ======== ============================= Support :math:`x \in (0, \infty)` - Mean :math:`dfrac{h}{4} if :math:`z=0`, :math:`\dfrac{tanh(z/2)h}{2z}` otherwise. + Mean :math:`\dfrac{h}{4}` if :math:`z=0`, :math:`\dfrac{tanh(z/2)h}{2z}` otherwise. Variance :math:`0.041666688h` if :math:`z=0`, :math:`\dfrac{h(sinh(z) - z)(1 - tanh^2(z/2))}{4z^3}` otherwise. ======== ============================= From 8129cbd588ce53da0b0747c417d1700e1967ba36 Mon Sep 17 00:00:00 2001 From: amyoshino Date: Tue, 20 Jun 2023 14:51:27 -0400 Subject: [PATCH 4/7] add icdf functions to halfcauchy and half normal --- pymc/distributions/continuous.py | 18 ++++++++++++++++++ tests/distributions/test_continuous.py | 8 ++++++++ 2 files changed, 26 insertions(+) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index 11ce4ae180..e6f5204284 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -856,6 +856,15 @@ def logcdf(value, loc, sigma): msg="sigma > 0", ) + def icdf(value, loc, sigma): + res = Normal.icdf((value + 1.0) / 2.0, loc, sigma) + res = check_icdf_value(res, value) + return check_icdf_parameters( + res, + sigma > 0, + msg="sigma > 0", + ) + class WaldRV(RandomVariable): name = "wald" @@ -2131,6 +2140,15 @@ def logcdf(value, loc, beta): msg="beta > 0", ) + def icdf(value, loc, beta): + res = loc + beta * pt.tan(np.pi * (value) / 2.0) + res = check_icdf_value(res, value) + return check_parameters( + res, + beta > 0, + msg="beta > 0", + ) + class Gamma(PositiveContinuous): r""" diff --git a/tests/distributions/test_continuous.py b/tests/distributions/test_continuous.py index 3c5dd5fe4c..dccf6afeff 100644 --- a/tests/distributions/test_continuous.py +++ b/tests/distributions/test_continuous.py @@ -299,6 +299,11 @@ def test_half_normal(self): {"sigma": Rplus}, lambda value, sigma: st.halfnorm.logcdf(value, scale=sigma), ) + check_icdf( + pm.HalfNormal, + {"sigma": Rplus}, + lambda q, sigma: st.halfnorm.ppf(q, scale=sigma), + ) def test_chisquared_logp(self): check_logp( @@ -582,6 +587,9 @@ def test_half_cauchy(self): {"beta": Rplusbig}, lambda value, beta: st.halfcauchy.logcdf(value, scale=beta), ) + check_icdf( + pm.HalfCauchy, {"beta": Rplusbig}, lambda q, beta: st.halfcauchy.ppf(q, scale=beta) + ) def test_gamma_logp(self): check_logp( From 63fd6bac7956f1f59535d27c392b757ac0432d87 Mon Sep 17 00:00:00 2001 From: amyoshino Date: Sat, 24 Jun 2023 07:39:40 -0400 Subject: [PATCH 5/7] obtaining icdf from full distribution --- pymc/distributions/continuous.py | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index e6f5204284..fc00e0f021 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -57,6 +57,7 @@ from pytensor.tensor.var import TensorConstant from pymc.logprob.abstract import _logcdf_helper, _logprob_helper +from pymc.logprob.basic import icdf try: from polyagamma import polyagamma_cdf, polyagamma_pdf, random_polyagamma @@ -857,13 +858,9 @@ def logcdf(value, loc, sigma): ) def icdf(value, loc, sigma): - res = Normal.icdf((value + 1.0) / 2.0, loc, sigma) + res = icdf(Normal.dist(loc, sigma), (value + 1.0) / 2.0) res = check_icdf_value(res, value) - return check_icdf_parameters( - res, - sigma > 0, - msg="sigma > 0", - ) + return res class WaldRV(RandomVariable): @@ -1731,13 +1728,9 @@ def logcdf(value, mu, sigma): ) def icdf(value, mu, sigma): - res = pt.exp(Normal.icdf(value, mu, sigma)) + res = pt.exp(icdf(Normal.dist(mu, sigma), value)) res = check_icdf_value(res, value) - return check_icdf_parameters( - res, - sigma > 0, - msg="sigma > 0", - ) + return res Lognormal = LogNormal From 08f6974700c850507a9cfdb76e1bff8d3497f92c Mon Sep 17 00:00:00 2001 From: amyoshino Date: Sat, 24 Jun 2023 20:44:05 -0400 Subject: [PATCH 6/7] retrigger tests From cdc104c4de900f86be3782e37606784ae7b9cd0a Mon Sep 17 00:00:00 2001 From: amyoshino Date: Sun, 25 Jun 2023 14:08:52 -0400 Subject: [PATCH 7/7] fixing lognormal --- pymc/distributions/continuous.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py index fc00e0f021..fa4c056f93 100644 --- a/pymc/distributions/continuous.py +++ b/pymc/distributions/continuous.py @@ -1729,7 +1729,6 @@ def logcdf(value, mu, sigma): def icdf(value, mu, sigma): res = pt.exp(icdf(Normal.dist(mu, sigma), value)) - res = check_icdf_value(res, value) return res