Add `loglogistic`, `logitexp`, `log1mlogistic` and `logit1mexp` #82

andrewjradcliffe · 2024-05-03T19:54:17Z

This takes advantage of LogExpFunctions's accurate implementations of log1pexp and log1mexp, combined with negation in the log-odds domain to provide more accurate and less expensive implementations of the function compositions.

As recommended here, I am submitting these.

Precedent for loglogistic and log1mlogistic is their existence in Stan; it is natural to include their inverses.

This takes advantage of `LogExpFunctions`'s accurate implementations of `log1pexp` and `log1mexp`, combined with negation in the log-odds domain to provide more accurate and less expensive implementations of the function compositions.

src/basicfuns.jl

devmotion · 2024-05-03T20:41:05Z

src/basicfuns.jl

+
+Its inverse is the [`loglogistic`](@ref) function.
+"""
+logitexp(x::Real) = x - log1mexp(x)


Should this be done more carefully to avoid loss of precision when both terms are approximately equal?

You're right -- this can be done more carefully using logexpm1.
However, as illustrated in the figures, logexpm1 itself suffers a loss of precision around -log(2).

Absolute error is ok.

But units in the last place error is a problem around -log(2)

As this can only be resolved by improving the ulp error of logexpm1 around log(2), I see no reason for this to gate merging of this PR.

I think we could check if the argument is in a neighborhood of -log(2), and then just use a first-order approximation. My quick calculations suggest that the local derivative is 2 but someone should check this.

src/basicfuns.jl

…expm1`

i.e. properly handle `Integer` and `Rational`

This enables the desired interface of `f(x::Real)`

devmotion · 2024-05-04T22:21:50Z

src/basicfuns.jl

+loglogistic(x::Real) = _loglogistic(x)
+_loglogistic(x::Real) = -log1pexp(-x)
+_loglogistic(x::Union{Integer, Rational}) = _loglogistic(float(x))


I'm still not sure whether unsigned integers are worth any consideration given that you will run into problems with them in many (most?) calculations and definitions in JuliaStats, JuliaMath, etc. typically do not handle them in a special way - but if we care then it would be easier to just define

Suggested change

loglogistic(x::Real) = _loglogistic(x)

_loglogistic(x::Real) = -log1pexp(-x)

_loglogistic(x::Union{Integer, Rational}) = _loglogistic(float(x))

loglogistic(x::Real) = -log1pexp(-float(x))

This would cover also cases such as e.g. unitful integer numbers that are not covered by Union{Integer, Rational} and for floating points such as x::Float64 or x::Float32 it would be equivalent to -log1pexp(-x).

I opted for the float solution, as this is simplest and, furthermore, replicates the handling of integer and rational arguments in Julia Base.

I'm still not sure whether unsigned integers are worth any consideration

SpecialFunctions.jl (another package I contributed to) converts most arguments via float, hence, it handles both the unsigned integer cases as well as the signed integer (i.e. typemin(T<:Signed)) edge case.

Use of float to convert integers and rationals to a floating point type for use in special functions is common throughout Base. Unless integers and rationals are converted, then a function interface defined as x::Real cannot be claimed. On the other hand, if the intention is to guarantee reasonable results only for floating point numbers, then this should clearly be stated in the package description. However, I am of the opinion that this package should follow the behavior of Base and SpecialFunctions, whereby integers and rationals are converted to floating point numbers internally.

Unless conversion is enforced, this leads to some obviously incorrect results, e.g., logistic(typemin(Int)) == 1.0, but, clearly, should be 0.0.

src/basicfuns.jl

devmotion · 2024-05-04T22:24:11Z

ext/LogExpFunctionsChainRulesCoreExt.jl

+ChainRulesCore.@scalar_rule(loglogistic(x::Real), (logistic(-x),))
+ChainRulesCore.@scalar_rule(log1mlogistic(x::Real), (-logistic(x),))
+ChainRulesCore.@scalar_rule(logitexp(x::Real), (inv(1 - exp(x)),))
+ChainRulesCore.@scalar_rule(logit1mexp(x::Real), (-inv(1 - exp(x)),))
+


IMO these are not needed - AD handles the functions already fine since rules for log1pexp etc. are defined.

Suggested change

ChainRulesCore.@scalar_rule(loglogistic(x::Real), (logistic(-x),))

ChainRulesCore.@scalar_rule(log1mlogistic(x::Real), (-logistic(x),))

ChainRulesCore.@scalar_rule(logitexp(x::Real), (inv(1 - exp(x)),))

ChainRulesCore.@scalar_rule(logit1mexp(x::Real), (-inv(1 - exp(x)),))

As directed, these definitions have been removed.

Looking through the issues, #36 seems to indicate that these functions should support ChainRulesCore forward and reverse rules.

No, that's just not correct. These new functions are supported by ChainRulesCore-based AD systems just fine without any additional definitions. See https://juliadiff.org/ChainRulesCore.jl/stable/rule_author/which_functions_need_rules.html for some additional information on this matter.

andrewjradcliffe · 2024-05-05T21:24:21Z

As

all tests (included the added tests) are passing
documentation has been updated
the extension interfaces modified to support these 4 functions
and tests of the package extensions have been added

please let me know of any additional changes which must be made for this PR to be merged.

tpapp

LGTM, as far as I can see all requests have been addressed. Thanks for the fine work!

andrewjradcliffe · 2024-05-07T02:02:41Z

May we clarify: include or exclude ChainRulesCore forward/reverse rules?

Last item to adjust, then can merge.

andrewjradcliffe · 2024-05-08T21:28:43Z

It would be best to not let this bit rot. Moreover, it is easiest to make changes when it is in recent memory.

@devmotion could you clarify with @tpapp whether to include or exclude ChainRulesCore forward/reverse rules?

Once clarified, I can either revert the last commit or leave it in place. Should I then assume that there would be no further obstacles to merging?

tpapp · 2024-05-09T08:19:40Z

Feel free to remove them. Thanks for the ping!

This reverts commit 1b0f914.

andrewjradcliffe · 2024-05-17T18:04:28Z

Reverted ChainRulesCore support. Should be ready for merge -- @devmotion and @tpapp

devmotion · 2024-05-17T19:27:19Z

test/basicfuns.jl

+    lim2 = T === Float16 ? -10.0 : -37.0
+    xs = T[Inf, -Inf, 0.0, lim1, lim2]
+    for x in xs
+        @test loglogistic(x) == log(logistic(x))


If we test multiple types T, I think it would be good to also check that - as desired - the results are of type T.

Makes sense. Comprehensive return type tests have been added for the 4 functions included in this PR.

Anything else, @devmotion?

andrewjradcliffe added 4 commits May 3, 2024 12:49

Register inverses; why not?

218d508

Update documentation

f89c9c5

Fix failing test: be slightly less stringent

33e6677

devmotion reviewed May 3, 2024

View reviewed changes

Re-style documentation

51b9a5b

devmotion reviewed May 3, 2024

View reviewed changes

src/basicfuns.jl Outdated Show resolved Hide resolved

andrewjradcliffe added 14 commits May 3, 2024 15:14

Use alternate design for loglogistic

9f24823

Simplify dispatch cases

a7f9cdd

Add ChangesOfVariables support

223327d

Fix duplication in docstring

27322fc

Be as careful as possible in the evaluation; we get ulp error of `log…

716b886

…expm1`

Apply the same conversion technique as loglogistic

3da3f8a

i.e. properly handle `Integer` and `Rational`

Eliminate 3 lines of code

8ee3db1

Hide the proper handling of integers and rationals.

c9ab3f3

This enables the desired interface of `f(x::Real)`

Add tests of inverses

ac17f3c

Add ChainRules extension and respective tests

16d6c9c

Add ChangesOfVariables tests

0d65559

Fix typo

5d84266

Update documentation

86ac5f5

Add tests to ensure correctness for subtypes of Real in Base

eb1633a

devmotion reviewed May 4, 2024

View reviewed changes

andrewjradcliffe added 5 commits May 5, 2024 13:09

Simplify definitions

46c0e60

Remove ChainRules definitions

18a58ef

Also remove ChainRules tests

519a6a4

Match expected form of ChangesOfVariables tests

716b671

Fix oversight in function type dispatches

58133ad

Add ChainRulesCore support back

1b0f914

tpapp approved these changes May 6, 2024

View reviewed changes

Revert "Add ChainRulesCore support back"

2aa9f44

This reverts commit 1b0f914.

devmotion reviewed May 17, 2024

View reviewed changes

andrewjradcliffe added 3 commits May 24, 2024 09:51

Add extensive return type tests

177cebb

Fix minor oversight from running tests locally

0d60abd

Add tests for NaN

15cb4ef

devmotion approved these changes May 31, 2024

View reviewed changes

devmotion merged commit 6fa1f71 into JuliaStats:master May 31, 2024
4 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add `loglogistic`, `logitexp`, `log1mlogistic` and `logit1mexp` #82

Add `loglogistic`, `logitexp`, `log1mlogistic` and `logit1mexp` #82

andrewjradcliffe commented May 3, 2024 •

edited

Loading

devmotion May 3, 2024

andrewjradcliffe May 4, 2024

andrewjradcliffe May 5, 2024

tpapp May 6, 2024

devmotion May 4, 2024

andrewjradcliffe May 5, 2024

devmotion May 4, 2024

andrewjradcliffe May 5, 2024

andrewjradcliffe May 6, 2024

devmotion May 6, 2024

andrewjradcliffe commented May 5, 2024 •

edited

Loading

tpapp left a comment

andrewjradcliffe commented May 7, 2024

andrewjradcliffe commented May 8, 2024

tpapp commented May 9, 2024

andrewjradcliffe commented May 17, 2024

devmotion May 17, 2024

andrewjradcliffe May 24, 2024

andrewjradcliffe May 28, 2024

Add loglogistic, logitexp, log1mlogistic and logit1mexp #82

Add loglogistic, logitexp, log1mlogistic and logit1mexp #82

Conversation

andrewjradcliffe commented May 3, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andrewjradcliffe commented May 5, 2024 • edited Loading

tpapp left a comment

Choose a reason for hiding this comment

andrewjradcliffe commented May 7, 2024

andrewjradcliffe commented May 8, 2024

tpapp commented May 9, 2024

andrewjradcliffe commented May 17, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Add `loglogistic`, `logitexp`, `log1mlogistic` and `logit1mexp` #82

Add `loglogistic`, `logitexp`, `log1mlogistic` and `logit1mexp` #82

andrewjradcliffe commented May 3, 2024 •

edited

Loading

andrewjradcliffe commented May 5, 2024 •

edited

Loading