Reapply #113, update compat entries, and prepare StatsFuns 1.0.0

.github/workflows/ci.yml

@@ -26,8 +26,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
-          - '1'
+          - '1.3' # oldest supported version
+          - '1.6' # LTS
+          - '1' # latest release
         os:
           - ubuntu-latest
           - macos-latest

Project.toml

@@ -1,13 +1,12 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.18"
+version = "1.0.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Rmath = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
@@ -16,10 +15,9 @@ ChainRulesCore = "1"
 InverseFunctions = "0.1"
 IrrationalConstants = "0.1"
 LogExpFunctions = "0.3.2"
-Reexport = "1"
-Rmath = "0.4, 0.5, 0.6, 0.7"
-SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2"
-julia = "1"
+Rmath = "0.6.1, 0.7"
+SpecialFunctions = "1.8.4, 2.1.4"
+julia = "1.3"
 
 [extras]
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"

README.md

@@ -5,60 +5,13 @@ Mathematical functions related to statistics.
 [![CI](https://github.com/JuliaStats/StatsFuns.jl/actions/workflows/ci.yml/badge.svg?branch=master)](https://github.com/JuliaStats/StatsFuns.jl/actions/workflows/ci.yml?query=branch%3Amaster)
 [![codecov](https://codecov.io/gh/JuliaStats/StatsFuns.jl/branch/master/graph/badge.svg?token=CV9hGTr6cW)](https://codecov.io/gh/JuliaStats/StatsFuns.jl)
 
-This package provides a collection of mathematical constants and numerical functions for statistical computing.
+This package provides a collection of numerical functions for statistical computing.
 
-## Constants
-
-```julia
-twoπ,       # 2π
-fourπ,      # 4π
-halfπ,      # π / 2
-quartπ,     # π / 4
-invπ,       # 1 / π
-twoinvπ,    # 2 / π
-fourinvπ,   # 4 / π
-inv2π,      # 1 / (2π)
-inv4π,      # 1 / (4π)
-sqrt2,      # √2
-sqrt3,      # √3
-sqrtπ,      # √π
-sqrt2π,     # √2π
-sqrt4π,     # √4π
-sqrthalfπ,  # √(π / 2)
-invsqrt2,   # 1 / √2
-invsqrt2π,  # 1 / √2π
-loghalf,    # log(1 / 2)
-logtwo,     # log(2)
-logπ,       # log(π)
-log2π,      # log(2π)
-log4π,      # log(4π)
-```
+Other mathematical constants and numerical functions are provided by the packages [IrrationalConstants.jl](https://github.com/JuliaMath/IrrationalConstants.jl) and [LogExpFunctions.jl](https://github.com/JuliaStats/LogExpFunctions.jl).
 
 ## Basic Functions
 
 ```julia
-# basicfuns
-xlogx,          # x * log(x), or 0 when x is zero
-xlogy,          # x * log(y), or 0 when x is zero
-xlog1py,        # x * log(1 + y) for x > 0, or 0 when x == 0
-logistic,       # 1 / (1 + exp(-x))
-logit,          # log(x / (1 - x))
-log1psq,        # log(1 + x^2)
-log1pexp,       # log(1 + exp(x))
-log1mexp,       # log(1 - exp(x))
-log2mexp,       # log(2 - exp(x))
-logexpm1,       # log(exp(x) - 1)
-softplus,       # alias of log1pexp
-invsoftplus,    # alias of logexpm1
-log1pmx,        # log(1 + x) - x
-logmxp1,        # log(x) - x + 1
-logaddexp,      # log(exp(x) + exp(y))
-logsubexp,      # log(abs(exp(x) - exp(y)))
-logsumexp,      # log(sum(exp(x)))
-softmax,        # exp(x_i) / sum(exp(x)), for i
-softmax!,       # inplace softmax
-
-# misc
 logmvgamma,     # logarithm of multivariate gamma function
 lstirling_asym
 ```

src/StatsFuns.jl

@@ -1,59 +1,12 @@
-__precompile__(true)
-
 module StatsFuns
 
 using Base: Math.@horner
-using Reexport
+using IrrationalConstants
+using LogExpFunctions
 using SpecialFunctions
 import ChainRulesCore
 import InverseFunctions
 
-# reexports
-@reexport using IrrationalConstants:
-    twoπ,       # 2π
-    fourπ,      # 4π
-    halfπ,      # π / 2
-    quartπ,     # π / 4
-    invπ,       # 1 / π
-    twoinvπ,    # 2 / π
-    fourinvπ,   # 4 / π
-    inv2π,      # 1 / (2π)
-    inv4π,      # 1 / (4π)
-    sqrt2,      # √2
-    sqrt3,      # √3
-    sqrtπ,      # √π
-    sqrt2π,     # √2π
-    sqrt4π,     # √4π
-    sqrthalfπ,  # √(π / 2)
-    invsqrt2,   # 1 / √2
-    invsqrt2π,  # 1 / √2π
-    loghalf,    # log(1 / 2)
-    logtwo,     # log(2)
-    logπ,       # log(π)
-    log2π,      # log(2π)
-    log4π       # log(4π)
-
-@reexport using LogExpFunctions:
-    xlogx,          # x * log(x) for x > 0, or 0 when x == 0
-    xlogy,          # x * log(y) for x > 0, or 0 when x == 0
-    xlog1py,        # x * log(1 + y) for x > 0, or 0 when x == 0
-    logistic,       # 1 / (1 + exp(-x))
-    logit,          # log(x / (1 - x))
-    log1psq,        # log(1 + x^2)
-    log1pexp,       # log(1 + exp(x))
-    log1mexp,       # log(1 - exp(x))
-    log2mexp,       # log(2 - exp(x))
-    logexpm1,       # log(exp(x) - 1)
-    softplus,       # alias of log1pexp
-    invsoftplus,    # alias of logexpm1
-    log1pmx,        # log(1 + x) - x
-    logmxp1,        # log(x) - x + 1
-    logaddexp,      # log(exp(x) + exp(y))
-    logsubexp,      # log(abs(e^x - e^y))
-    logsumexp,      # log(sum(exp(x)))
-    softmax,        # exp(x_i) / sum(exp(x)), for i
-    softmax!        # inplace softmax
-
 export
     # distrs/beta
     betapdf,            # pdf of beta distribution

test/chainrules.jl

@@ -6,7 +6,7 @@ using Random
 @testset "chainrules" begin
     x = exp(randn())
     y = exp(randn())
-    z = logistic(randn())
+    z = StatsFuns.logistic(randn())
     test_frule(betalogpdf, x, y, z)
     test_rrule(betalogpdf, x, y, z)
 
@@ -33,7 +33,7 @@ using Random
     test_rrule(tdistlogpdf, x, y)
 
     x = rand(1:100)
-    y = logistic(randn())
+    y = StatsFuns.logistic(randn())
     z = rand(1:x)
     test_frule(binomlogpdf, x, y, z)
     test_rrule(binomlogpdf, x, y, z)

test/misc.jl

@@ -12,15 +12,15 @@ using SpecialFunctions, StatsFuns
         #  Γ₁(a) = Γ(a), Γ₂(a) = √π Γ(a) Γ(a - 0.5), etc
         a = rand(eltya) + 1 # add one since loggamma is only define for positive arguments
         @test logmvgamma(1, a) ≈ loggamma(a)
-        @test logmvgamma(2, a) ≈    eltya(0.5logπ) + loggamma(a) + loggamma(a - eltya(0.5))
-        @test logmvgamma(3, a) ≈ eltya((3/2)*logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a))
+        @test logmvgamma(2, a) ≈    eltya(0.5*StatsFuns.logπ) + loggamma(a) + loggamma(a - eltya(0.5))
+        @test logmvgamma(3, a) ≈ eltya((3/2)*StatsFuns.logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a))
     end
 
     @testset "consistent with itself" for eltya in (Float32, Float64)
         #  Γᵢ(a) = (π^{i-1/2}) Γ(a) Γᵢ₋₁(a - 0.5)
         for p in 1:50
             a = rand(eltya) + p # add p since loggamma is only define for positive arguments
-            @test logmvgamma(p, a) ≈ eltya((p/2-1/2)*logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5))
+            @test logmvgamma(p, a) ≈ eltya((p/2-1/2)*StatsFuns.logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5))
         end
     end
 end

test/rmath.jl

@@ -140,26 +140,24 @@ end
         ((10, 2), 0//1:1//100:1//1),
     ])
 
-    # We test the following extreme parameters separately since
-    # a slightly larger tolerance is needed.
-    #
-    # For `betapdf(1000, 2, 0.58)`:
-    # StatsFuns:   1.9419987107407202e-231
-    # Rmath:       1.941998710740941e-231
-    # Mathematica: 1.941998710742487e-231
-    # For `betapdf(1000, 2, 0.68)`:
-    # StatsFuns:   1.5205049885199752e-162
-    # Rmath:       1.5205049885200616e-162
-    # Mathematica: 1.520504988521358e-162
-    #
-    # With older SpecialFunctions + Julia versions a slightly larger tolerance is needed
-    @testset "Beta(1000, 2)" begin
-        if VERSION < v"1.3"
-            rmathcomp("beta", (1000, 2), 0.0:0.01:1.0, 1e-12)
-        else
-            rmathcomp("beta", (1000, 2), setdiff(0.0:0.01:1.0, (0.58, 0.68)))
-            rmathcomp("beta", (1000, 2), [0.58, 0.68], 1e-12)
-        end
+    # It is not possible to maintain a rtol of 1e-14 for the cdf when there is such a large difference
+    # in the magnitude of the parameters. At least not with the current parameters. Furthermore, it
+    # seems that Rmath is actually less accurate than we are but since we are comparing against Rmath
+    # we have to use rtol=1e-12 although we are probably only off by around 1e-13.
+    @testset "param: (1000, 2)" begin
+        rmathcomp(
+            "beta",
+            (1000, 2),
+            # We have to drop the 0.48 value since the R quantile function fails while we succeed.
+            # It's tested separate below.
+            setdiff(collect(0.0:0.01:1.0), 0.48),
+            1e-12)
+        # Test p=0.48 separately since R fails. (It's pretty slow, though, caused by the cdf being 9.0797754e-317)
+        @test betainvcdf(1000, 2, betacdf(1000, 2, 0.48)) ≈ 0.48
+    end
+    @testset "$(f)(0, 0, $x) should error" for f in (betacdf, betaccdf, betalogcdf, betalogccdf),
+        x in (0.0, 0.5, 1.0)
+        @test_throws DomainError f(0.0, 0.0, x)
     end
 
     rmathcomp_tests("binom", [