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

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"

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)