From c3fcc96dd20e0378648a671963fff4c52e298098 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 29 Dec 2020 18:09:00 +0100 Subject: [PATCH 1/5] Fix cutoff value in `log1mexp` and redundant reimplementation in `Exponential.logcdf()` --- pymc3/distributions/continuous.py | 12 +++--------- pymc3/math.py | 4 ++-- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 08102d2022..a10c158f45 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -45,7 +45,7 @@ ) from pymc3.distributions.distribution import Continuous, draw_values, generate_samples from pymc3.distributions.special import log_i0 -from pymc3.math import invlogit, logdiffexp, logit +from pymc3.math import invlogit, log1mexp, logdiffexp, logit from pymc3.theanof import floatX __all__ = [ @@ -1513,12 +1513,6 @@ def logcdf(self, value): Compute the log of cumulative distribution function for the Exponential distribution at the specified value. - References - ---------- - .. [Machler2012] Martin Mächler (2012). - "Accurately computing :math:`\log(1-\exp(-\mid a \mid))` Assessed by the Rmpfr - package" - Parameters ---------- value: numeric @@ -1533,9 +1527,9 @@ def logcdf(self, value): lam = self.lam a = lam * value return tt.switch( - tt.le(value, 0.0), + tt.or_(tt.le(value, 0.0), tt.le(lam, 0)), -np.inf, - tt.switch(tt.le(a, tt.log(2.0)), tt.log(-tt.expm1(-a)), tt.log1p(-tt.exp(-a))), + log1mexp(a), ) diff --git a/pymc3/math.py b/pymc3/math.py index 17f286ceab..9dddd67a0b 100644 --- a/pymc3/math.py +++ b/pymc3/math.py @@ -226,7 +226,7 @@ def log1mexp(x): For details, see https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf """ - return tt.switch(tt.lt(x, 0.683), tt.log(-tt.expm1(-x)), tt.log1p(-tt.exp(-x))) + return tt.switch(tt.lt(x, 0.693147), tt.log(-tt.expm1(-x)), tt.log1p(-tt.exp(-x))) def log1mexp_numpy(x): @@ -235,7 +235,7 @@ def log1mexp_numpy(x): For details, see https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf """ - return np.where(x < 0.683, np.log(-np.expm1(-x)), np.log1p(-np.exp(-x))) + return np.where(x < 0.693147, np.log(-np.expm1(-x)), np.log1p(-np.exp(-x))) def flatten_list(tensors): From 7ca4f31d8fd9a389c8017c95de0c80493a747cc2 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 29 Dec 2020 19:15:58 +0100 Subject: [PATCH 2/5] Get more digits :) --- pymc3/math.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc3/math.py b/pymc3/math.py index 9dddd67a0b..c2d1ae1555 100644 --- a/pymc3/math.py +++ b/pymc3/math.py @@ -226,7 +226,7 @@ def log1mexp(x): For details, see https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf """ - return tt.switch(tt.lt(x, 0.693147), tt.log(-tt.expm1(-x)), tt.log1p(-tt.exp(-x))) + return tt.switch(tt.lt(x, 0.6931471805599453), tt.log(-tt.expm1(-x)), tt.log1p(-tt.exp(-x))) def log1mexp_numpy(x): @@ -235,7 +235,7 @@ def log1mexp_numpy(x): For details, see https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf """ - return np.where(x < 0.693147, np.log(-np.expm1(-x)), np.log1p(-np.exp(-x))) + return np.where(x < 0.6931471805599453, np.log(-np.expm1(-x)), np.log1p(-np.exp(-x))) def flatten_list(tensors): From 2af5b03891ea260bfb2d568e82b8a1333b1b1ab3 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Thu, 31 Dec 2020 13:25:20 +0100 Subject: [PATCH 3/5] Remove redundant reimplementation in Weibull --- pymc3/distributions/continuous.py | 11 +++-------- pymc3/math.py | 9 ++++++++- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index a10c158f45..cb156367ef 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -2800,12 +2800,6 @@ def logcdf(self, value): Compute the log of the cumulative distribution function for Weibull distribution at the specified value. - References - ---------- - .. [Machler2012] Martin Mächler (2012). - "Accurately computing `\log(1-\exp(- \mid a \mid))` Assessed by the Rmpfr - package" - Parameters ---------- value: numeric @@ -2822,7 +2816,7 @@ def logcdf(self, value): return tt.switch( tt.le(value, 0.0), -np.inf, - tt.switch(tt.le(a, tt.log(2.0)), tt.log(-tt.expm1(-a)), tt.log1p(-tt.exp(-a))), + log1mexp(a), ) @@ -3902,7 +3896,7 @@ def logcdf(self, value): References ---------- - .. [Machler2012] Martin Mächler (2012). + .. [Machler2012] c. "Accurately computing :math: `\log(1-\exp(- \mid a \mid<))` Assessed by the Rmpfr package" @@ -3916,6 +3910,7 @@ def logcdf(self, value): ------- TensorVariable """ + # TODO: Check possible redundant (or improved) reimplementation of log1pexp mu = self.mu s = self.s a = -(value - mu) / s diff --git a/pymc3/math.py b/pymc3/math.py index c2d1ae1555..17a7a7f67d 100644 --- a/pymc3/math.py +++ b/pymc3/math.py @@ -219,12 +219,19 @@ def log1pexp(x): def log1mexp(x): - """Return log(1 - exp(-x)). + r"""Return log(1 - exp(-x)). This function is numerically more stable than the naive approach. For details, see https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf + + References + ---------- + .. [Machler2012] Martin Mächler (2012). + "Accurately computing `\log(1-\exp(- \mid a \mid))` Assessed by the Rmpfr + package" + """ return tt.switch(tt.lt(x, 0.6931471805599453), tt.log(-tt.expm1(-x)), tt.log1p(-tt.exp(-x))) From fea4a87088b02b5eaa3d5bad3a3e77fc3ee0c822 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Sat, 2 Jan 2021 12:03:26 +0100 Subject: [PATCH 4/5] Tiny rewrite to use `|` instead of `tt.or_`. Remove TODO comments. --- pymc3/distributions/continuous.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index cb156367ef..2d023dbfbf 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -1527,7 +1527,7 @@ def logcdf(self, value): lam = self.lam a = lam * value return tt.switch( - tt.or_(tt.le(value, 0.0), tt.le(lam, 0)), + tt.le(value, 0.0) | tt.le(lam, 0), -np.inf, log1mexp(a), ) @@ -3896,7 +3896,7 @@ def logcdf(self, value): References ---------- - .. [Machler2012] c. + .. [Machler2012] Martin Mächler (2012). "Accurately computing :math: `\log(1-\exp(- \mid a \mid<))` Assessed by the Rmpfr package" @@ -3910,7 +3910,6 @@ def logcdf(self, value): ------- TensorVariable """ - # TODO: Check possible redundant (or improved) reimplementation of log1pexp mu = self.mu s = self.s a = -(value - mu) / s From d7f82a85dce9e0874f73dd6d48b096d5f8792313 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Sat, 2 Jan 2021 13:01:51 +0100 Subject: [PATCH 5/5] retrigger checks