Update log1mexp and remove redundant local reimplementations in the library

[ricardoV94]

The cutoff in log1mexp (and its numpy counterpart) is supposed to be log(2) = 0.693. There was probably a small typo when writing the function as the cutoff was set to 0.683. Pinging @aseyboldt

Original source of this algorithm: https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
Example R implementation: https://github.com/cran/Rmpfr/blob/7f15be8d124a7517719012107c132b5d828284f9/inst/doc/log1mexp-note.R#L96
TensorFlow implementation: https://github.com/tensorflow/probability/blob/v0.11.1/tensorflow_probability/python/math/generic.py#L495-L519
Haskell implementation: https://gitlab.haskell.org/ghc/ghc/-/blob/master/libraries/base/GHC/Float.hs#L159-161

I also removed a redundant re-implementation of this function in the logcdf method of the Exponential and Weibull distribution.

---

I would further suggest that the log1mexp method be changed to expect negative values (such as logarithms of probabilities) instead of positive values. This is what the Haskell implementation does directly and what the TensorFlow does less directly by first taking the absolute of the input and then negating it. The function would change from log(1 - exp(-x)) to log(1 - exp(x)) (as in Haskell). This is equivalent to log(1 - exp(-|x|)) (as in TFP), but I like the formula with the absolute less because it seems to just hide mistakes such as asking log(1 - exp(log(1.1)). I don't see why one would ever want values larger than 1 after exponentiation to be negated before exponentiation.

I have actually implemented this change and it did not seem to have large ramifications. There seem to be only 2 places (3 after #4387) in the library where the function is used directly, and the only thing needed is to negate the input (e.g., the function logdiffexp needs to be changed from a + log1mexp(a - b) to a + log1mexp(b - a)). I also checked that current unittests would detect all the spotted sections if they were not corrected.

The least nice part is that it would probably require a UserWarning about the change when the function is called directly by the user. What do you think?

Update: This change was rejected to avoid backwards incompatibility

[kyleabeauchamp]

What about people using double precision? I bet this constant ought to have about twice as many digits...

[ricardoV94]

Pardon my ignorance, how many digits would that be :p

My ipython prints: np.log(2) = 0.6931471805599453
We could also put directly np.log(2) or is that too wasteful?

Edit: I think I misunderstood your question, but anyway adding more digits seems like a no-brainer.

[codecov]

Codecov Report

Merging #4394 (d7f82a8) into master (3cfee77) will increase coverage by 0.05%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master    #4394      +/-   ##
==========================================
+ Coverage   88.04%   88.09%   +0.05%     
==========================================
  Files          88       88              
  Lines       14482    14538      +56     
==========================================
+ Hits        12750    12807      +57     
+ Misses       1732     1731       -1     

Impacted Files	Coverage Δ
pymc3/distributions/continuous.py	94.37% <100.00%> (ø)
pymc3/math.py	68.64% <100.00%> (ø)
pymc3/distributions/discrete.py	95.00% <0.00%> (-0.85%)	⬇️
pymc3/sampling_jax.py	0.00% <0.00%> (ø)
pymc3/distributions/multivariate.py	83.64% <0.00%> (+0.72%)	⬆️

[AlexAndorra]

Changes look good, thanks @ricardoV94 👌

Regarding this:

I would further suggest that the log1mexp method be changed to expect negative values [...] The least nice part is that it would probably require a UserWarning about the change when the function is called directly by the user. What do you think?

This change is not implemented in the PR yet, right? Or I didn't see it 🤔
If we do make this change though, a UserWarning is indeed necessary

[ricardoV94]

This change is not implemented in the PR yet, right? Or I didn't see it 🤔

No it's not, I wanted to ask other's opinion first.

[twiecki]

@ricardoV94 Let's go with your preference.

[ricardoV94]

@ricardoV94 Let's go with your preference.

Cool. I will push the changes to this PR

[ricardoV94]

TODO:

• Remove newly found redundant reimplementation from Weibull.logcdf()
• Rewrite log1mexp from log(1 - exp(-x)) to log(1 - exp(x))
• Add UserWarning about change in behavior

About the last point, what is the best way to avoid UserWarnings from internal API calls to the new log1mexp? Having a WarningFilter before every call seems inelegant. Perhaps write the new function in _log1mexp and have the warning only in log1mexp followed by a call to the new function? This way API calls can directly call _log1mexp without having to worry about the warning.

What do you think?

[twiecki]

@ricardoV94 This fixes a bug which we don't need to warn about, and adds support for negative inputs which was previously not allowed, which we also don't need to warn about. Are my assumptions incorrect?

[ricardoV94]

@ricardoV94 This fixes a bug which we don't need to warn about, and adds support for negative inputs which was previously not allowed, which we also don't need to warn about. Are my assumptions incorrect?

Yes and no.

The PR fixes a small inefficiency in the original method due to using a slightly off value 0.683... instead of 0.693... It would also remove two redundant reimplementations that I found in the codebase. Those things indeed do not need any warnings.

In addition to this, I suggested changing the sign of the expected input (purely on the ground that is is more intuitive). Right now it expects positive numbers which are explicitly negated, and I suggested expecting negative values instead and not negating them anymore. This would break things retroactively: it would no longer work with positive inputs (now it does not work with negative inputs). The other option is to use the TFP approach of negating the absolute. I find it less appealing but it has the advantage of being backwards compatible. However, there is nothing fundamentally wrong with the current approach and the documentation is correct.

[twiecki]

@ricardoV94 Thanks for explaining, I did not understand the full scope before. In that case I think we should stick with the previous one in the interest of backwards compat.

[ricardoV94]

@ricardoV94 Thanks for explaining, I did not understand the full scope before. In that case I think we should stick with the previous one in the interest of backwards compat.

Ok that's fine. I will just remove the other duplicated implementation that I found and this should be ready to merge.

[ricardoV94]

Ready for review / merging :)

I think this one does not warrant a release-note. Do you agree?

[AlexAndorra]

Looking good, thanks @ricardoV94 ! No need for a release note indeed I think, as it's really internal maintenance. I'll merge once tests pass 😉

[twiecki]

Thanks @ricardoV94!