Fix cutoff value in log1mexp and redundant reimplementation in Exponential.logcdf()

Author: ricardoV94

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),
         )
 
 

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):