Add loglogistic, logitexp, log1mlogistic and logit1mexp

docs/src/index.md

@@ -4,7 +4,7 @@ Various special functions based on `log` and `exp` moved from [StatsFuns.jl](htt
 
 The original authors of these functions are the StatsFuns.jl contributors.
 
-LogExpFunctions supports [`InverseFunctions.inverse`](https://github.com/JuliaMath/InverseFunctions.jl) and [`ChangesOfVariables.test_with_logabsdet_jacobian`](https://github.com/JuliaMath/ChangesOfVariables.jl) for `log1mexp`, `log1pexp`, `log2mexp`, `logexpm1`, `logistic`, `logit`, and `logcosh` (no inverse).
+LogExpFunctions supports [`InverseFunctions.inverse`](https://github.com/JuliaMath/InverseFunctions.jl) and [`ChangesOfVariables.test_with_logabsdet_jacobian`](https://github.com/JuliaMath/ChangesOfVariables.jl) for `log1mexp`, `log1pexp`, `log2mexp`, `logexpm1`, `logistic`, `logit`, `loglogistic`, `logitexp`, `log1mlogistic`, `logit1mexp`, and `logcosh` (no inverse).
 
 ```@docs
 xlogx
@@ -31,4 +31,8 @@ softmax!
 softmax
 cloglog
 cexpexp
+loglogistic
+logitexp
+log1mlogistic
+logit1mexp
 ```

ext/LogExpFunctionsChainRulesCoreExt.jl

@@ -187,4 +187,9 @@ end
 ChainRulesCore.@scalar_rule(cloglog(x), (-inv((1 - x) * log1p(-x)),))
 ChainRulesCore.@scalar_rule(cexpexp(x), (-xexpx(-exp(x)),))
 
+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)),))
+
 end # module

ext/LogExpFunctionsChangesOfVariablesExt.jl

@@ -42,4 +42,24 @@ function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logcosh), x::Real)
     return y, log1mexp(a) - z
 end
 
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(loglogistic), x::Real)
+    y = loglogistic(x)
+    return y, y - x
+end
+
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(log1mlogistic), x::Real)
+    y = log1mlogistic(x)
+    return y, x + y
+end
+
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logitexp), x::Real)
+    y = logitexp(x)
+    return y, y - x
+end
+
+function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logit1mexp), x::Real)
+    y = logit1mexp(x)
+    return y, -y - x
+end
+
 end # module

ext/LogExpFunctionsInverseFunctionsExt.jl

@@ -16,4 +16,10 @@ InverseFunctions.inverse(::typeof(logistic)) = logit
 InverseFunctions.inverse(::typeof(cloglog)) = cexpexp
 InverseFunctions.inverse(::typeof(cexpexp)) = cloglog
 
+InverseFunctions.inverse(::typeof(loglogistic)) = logitexp
+InverseFunctions.inverse(::typeof(logitexp)) = loglogistic
+
+InverseFunctions.inverse(::typeof(log1mlogistic)) = logit1mexp
+InverseFunctions.inverse(::typeof(logit1mexp)) = log1mlogistic
+
 end # module

src/LogExpFunctions.jl

@@ -8,7 +8,8 @@ import LinearAlgebra
 
 export xlogx, xlogy, xlog1py, xexpx, xexpy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1,
     softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, logsumexp!, softmax,
-    softmax!, logcosh, logabssinh, cloglog, cexpexp
+    softmax!, logcosh, logabssinh, cloglog, cexpexp,
+    loglogistic, logitexp, log1mlogistic, logit1mexp
 
 # expm1(::Float16) is not defined in older Julia versions,
 # hence for better Float16 support we use an internal function instead

src/basicfuns.jl

@@ -441,3 +441,68 @@ $(SIGNATURES)
 Compute the complementary double exponential, `1 - exp(-exp(x))`.
 """
 cexpexp(x) = -_expm1(-exp(x))
+
+#=
+this uses the identity:
+
+log(logistic(x)) = -log(1 + exp(-x))
+=#
+"""
+$(SIGNATURES)
+
+Return `log(logistic(x))`, computed more carefully and with fewer calls
+than the naive composition of functions.
+
+Its inverse is the [`logitexp`](@ref) function.
+"""
+loglogistic(x::Real) = -log1pexp(-float(x))
+
+#=
+this uses the identity:
+
+logit(exp(x)) = log(exp(x) / (1 + exp(x))) = -log(exp(-x) - 1)
+=#
+"""
+$(SIGNATURES)
+
+Return `logit(exp(x))`, computed more carefully and with fewer calls than
+the naive composition of functions.
+
+Its inverse is the [`loglogistic`](@ref) function.
+"""
+logitexp(x::Real) = -logexpm1(-float(x))
+
+#=
+this uses the identity:
+
+log(logistic(-x)) = -log(1 + exp(x))
+
+that is, negation in the log-odds domain.
+=#
+
+"""
+$(SIGNATURES)
+
+Return `log(1 - logistic(x))`, computed more carefully and with fewer calls than
+the naive composition of functions.
+
+Its inverse is the [`logit1mexp`](@ref) function.
+"""
+log1mlogistic(x::Real) = -log1pexp(x)
+
+#=
+
+this uses the same identity:
+
+-logit(exp(x)) = logit(1 - exp(x)) = log((1 - exp(x)) / exp(x)) = log(exp(-x) - 1)
+=#
+
+"""
+$(SIGNATURES)
+
+Return `logit(1 - exp(x))`, computed more carefully and with fewer calls than
+the naive composition of functions.
+
+Its inverse is the [`log1mlogistic`](@ref) function.
+"""
+logit1mexp(x::Real) = logexpm1(-float(x))

test/basicfuns.jl

@@ -468,3 +468,93 @@ end
     @test cexpexp(-Inf) == 0.0
     @test cexpexp(0) == (ℯ - 1) / ℯ
 end
+
+@testset "loglogistic: $T" for T in (Float16, Float32, Float64)
+    lim1 = T === Float16 ? -14.0 : -50.0
+    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))
+    end
+
+    ϵ = eps(T)
+    xs = T[ϵ, 1.0, 18.0, 33.3, 50.0]
+    for x in xs
+        lhs = loglogistic(x)
+        rhs = log(logistic(x))
+        @test abs(lhs - rhs) < ϵ
+    end
+
+    # misc
+    @test loglogistic(T(Inf)) == -zero(T)
+    @test loglogistic(-T(Inf)) == -T(Inf)
+    @test loglogistic(-T(103.0)) == -T(103.0)
+    @test abs(loglogistic(T(35.0))) < 3eps(T)
+    @test abs(loglogistic(T(103.0))) < eps(T)
+    @test isfinite(loglogistic(-T(745.0)))
+    @test isfinite(loglogistic(T(50.0)))
+    @test isfinite(loglogistic(T(745.0)))
+end
+
+
+@testset "logitexp: $T" for T in (Float16, Float32, Float64)
+    ϵ = eps(T)
+    xs = T[ϵ, √ϵ, 0.2, 0.4, 0.8, 1.0 - √ϵ, 1.0 - ϵ]
+    neg_xs = -xs
+    for x in xs
+        @test abs(logitexp(loglogistic(x)) - x) <= ϵ
+    end
+    for x in neg_xs
+        @test abs(logitexp(loglogistic(x)) - x) <= 2ϵ
+    end
+    xs = T[-Inf, 0.0, Inf]
+    for x in xs
+        @test logitexp(loglogistic(x)) == x
+    end
+end
+
+@testset "log1mlogistic: $T" for T in (Float16, Float32, Float64)
+    @test log1mlogistic(T(Inf)) == -T(Inf)
+    @test log1mlogistic(-T(Inf)) == -zero(T)
+    @test log1mlogistic(-T(103.0)) < eps(T)
+    @test abs(log1mlogistic(T(35.0))) == T(35.0)
+    @test abs(log1mlogistic(T(103.0))) == T(103.0)
+    @test isfinite(log1mlogistic(-T(745.0)))
+    @test isfinite(log1mlogistic(T(50.0)))
+    @test isfinite(log1mlogistic(T(745.0)))
+end
+
+
+@testset "logit1mexp: $T" for T in (Float16, Float32, Float64)
+    ϵ = eps(T)
+    xs = T[ϵ, √ϵ, 0.2, 0.4, 0.8, 1.0 - √ϵ, 1.0 - ϵ]
+    neg_xs = -xs
+    for x in xs
+        @test abs(logit1mexp(log1mlogistic(x)) - x) <= 2ϵ
+    end
+    for x in neg_xs
+        @test abs(logit1mexp(log1mlogistic(x)) - x) <= ϵ
+    end
+    xs = T[-Inf, 0.0, Inf]
+    for x in xs
+        @test logit1mexp(log1mlogistic(x)) == x
+    end
+end
+
+@testset "correctness wrt Unsigned, Rational" begin
+    @test loglogistic(UInt64(5)) == loglogistic(5.0)
+    @test log1mlogistic(UInt64(5)) == log1mlogistic(5.0)
+
+    @test loglogistic(0x01//0x02) == loglogistic(0.5)
+    @test log1mlogistic(0x01//0x02) == log1mlogistic(0.5)
+end
+
+@testset "correctness wrt Integer edge case" begin
+    # If not handled, these will be zero
+    @test loglogistic(typemin(Int)) == loglogistic(float(typemin(Int)))
+    @test log1mlogistic(typemin(Int)) == log1mlogistic(float(typemin(Int)))
+
+    # If not handled, these would throw since negation at typemin is a round-trip
+    @test logitexp(typemin(Int)) == logitexp(float(typemin(Int)))
+    @test logit1mexp(typemin(Int)) == logit1mexp(float(typemin(Int)))
+end

test/chainrules.jl

@@ -125,4 +125,30 @@
     test_scalar(cloglog, rand())
 
     test_scalar(cexpexp, rand())
+
+    @testset "loglogistic, log1mlogistic" begin
+        # relationship to logistic suggests an analogous test regime
+        for x in (-821.4, -23.5, 12.3, 41.2)
+            test_frule(loglogistic, x)
+            test_rrule(loglogistic, x)
+        end
+        for x in (-123.2f0, -21.4f0, 8.3f0, 21.5f0)
+            test_frule(loglogistic, x; rtol=1f-3, atol=1f-3)
+            test_rrule(loglogistic, x; rtol=1f-3, atol=1f-3)
+        end
+
+        for x in (-821.4, -23.5, 12.3, 41.2)
+            test_frule(log1mlogistic, x)
+            test_rrule(log1mlogistic, x)
+        end
+        for x in (-123.2f0, -21.4f0, 8.3f0, 21.5f0)
+            test_frule(log1mlogistic, x; rtol=1f-3, atol=1f-3)
+            test_rrule(log1mlogistic, x; rtol=1f-3, atol=1f-3)
+        end
+    end
+
+    test_frule(logitexp, -x)
+    test_frule(logit1mexp, -x)
+    test_rrule(logitexp, -x)
+    test_rrule(logit1mexp, -x)
 end

test/inverse.jl

@@ -11,4 +11,10 @@
 
     InverseFunctions.test_inverse(cloglog, rand())
     InverseFunctions.test_inverse(cexpexp, rand())
+
+    InverseFunctions.test_inverse(loglogistic, randexp())
+    InverseFunctions.test_inverse(logitexp, -randexp())
+
+    InverseFunctions.test_inverse(log1mlogistic, randexp())
+    InverseFunctions.test_inverse(logit1mexp, -randexp())
 end

test/with_logabsdet_jacobian.jl

@@ -19,4 +19,10 @@
 
     ChangesOfVariables.test_with_logabsdet_jacobian(logcosh, x, derivative)
     ChangesOfVariables.test_with_logabsdet_jacobian(logcosh, -x, derivative)
+
+    ChangesOfVariables.test_with_logabsdet_jacobian(loglogistic, x, derivative)
+    ChangesOfVariables.test_with_logabsdet_jacobian(logitexp, -x, derivative)
+
+    ChangesOfVariables.test_with_logabsdet_jacobian(log1mlogistic, x, derivative)
+    ChangesOfVariables.test_with_logabsdet_jacobian(logit1mexp, -x, derivative)
 end