newton(): defer to Roots.jl implementation

Note that we really do need `Roots.jl` 2.2.0-or-newer,
because of JuliaMath/Roots.jl#445

This forces minimal Julia version to be 1.6

.github/workflows/CI.yml

@@ -23,7 +23,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
+          - '1.6'
           - '1'
           - pre
         os:
@@ -36,7 +36,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           # ARM64 on macos-latest is neither supported by older Julia versions nor setup-julia
-          arch: ${{ matrix.os == 'macos-latest' && matrix.version != '1.3' && 'aarch64' || 'x64' }}
+          arch: ${{ matrix.os == 'macos-latest' && matrix.version >= '1.8' && 'aarch64' || 'x64' }}
           show-versioninfo: true
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1

Project.toml

@@ -13,6 +13,7 @@ PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
@@ -48,6 +49,7 @@ PDMats = "0.10, 0.11"
 Printf = "<0.0.1, 1"
 QuadGK = "2"
 Random = "<0.0.1, 1"
+Roots = "2.2"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "1.2, 2"
 StableRNGs = "1"
@@ -57,7 +59,7 @@ StatsAPI = "1.6"
 StatsBase = "0.32, 0.33, 0.34"
 StatsFuns = "0.9.15, 1"
 Test = "<0.0.1, 1"
-julia = "1.3"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

src/quantilealgs.jl

@@ -1,3 +1,5 @@
+using Roots
+
 # Various algorithms for computing quantile
 
 function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx::T) where {T<:Real}
@@ -47,20 +49,8 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
-function newton_impl(Δ, xs::T=mode(d), xrtol::Real=1e-12) where {T}
-    x = xs - Δ(xs)
-    @assert typeof(x) === T
-    x0 = T(xs)
-    while !isapprox(x, x0, atol=0, rtol=xrtol)
-        x0 = x
-        x = x0 - Δ(x0)
-    end
-    return x
-end
-
 function newton((f,df), xs::T=mode(d), xrtol::Real=1e-12) where {T}
-    Δ(x) = f(x)/df(x)
-    return newton_impl(Δ, xs, xrtol)
+    return find_zero((f,df), xs, Roots.Newton(), xatol=0, xrtol=xrtol, atol=0, rtol=eps(float(T)), maxiters=typemax(Int))
 end
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)