diff --git a/Project.toml b/Project.toml
index 092dfae81..bb91ff1f6 100644
--- a/Project.toml
+++ b/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"
@@ -51,6 +52,7 @@ PDMats = "0.11.35"
 Printf = "<0.0.1, 1"
 QuadGK = "2"
 Random = "<0.0.1, 1"
+Roots = "2"
 SparseArrays = "<0.0.1, 1"
 SparseConnectivityTracer = "1"
 SpecialFunctions = "1.2, 2"
diff --git a/src/Distributions.jl b/src/Distributions.jl
index ed2a2cdde..d00daf121 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -4,6 +4,7 @@ using StatsBase, PDMats, StatsFuns, Statistics
 using StatsFuns: logtwo, invsqrt2, invsqrt2π
 
 import QuadGK: quadgk
+using Roots: find_zero, ITP, Newton
 import Base: size, length, convert, show, getindex, rand, vec, inv
 import Base: sum, maximum, minimum, extrema, +, -, *, ==
 import Base.Math: @horner
diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 8aa9f1b89..8976cab07 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -3,35 +3,20 @@
 function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx::T) where {T<:Real}
     rx < lx && throw(ArgumentError("empty bracketing interval [$lx, $rx]"))
 
-    # In some special cases, e.g. #1501, rx == lx`
-    # If the distribution is degenerate the check below can fail, hence we skip it
-    if rx == lx
-        # Returns `lx` of the same type as `(lx + rx) / 2`
-        # For specific types such as `Float64` it is more performant than `oftype((lx + rx) / 2, lx)`
-        return middle(lx)
-    end
-
     # base tolerance on types to support e.g. `Float32` (avoids an infinite loop)
     # ≈ 3.7e-11 for Float64
     # ≈ 2.4e-5 for Float32
     tol = cbrt(eps(float(T)))^2
-    # find quantile using bisect algorithm
-    cl = cdf(d, lx)
-    cr = cdf(d, rx)
-    cl <= p <= cr ||
-        throw(ArgumentError("[$lx, $rx] is not a valid bracketing interval for `quantile(d, $p)`"))
-    while rx - lx > tol
-        m = (lx + rx)/2
-        c = cdf(d, m)
-        if p > c
-            cl = c
-            lx = m
-        else
-            cr = c
-            rx = m
-        end
+
+    # In some special cases, e.g. #1501, rx == lx`
+    # If the distribution is degenerate the check below can fail, hence we skip it
+    if lx == rx || rx - lx <= tol * max(abs(lx), abs(rx))
+        # Returns `lx` of the same type as `(lx + rx) / 2`
+        # For specific types such as `Float64` it is more performant than `oftype((lx + rx) / 2, lx)`
+        return middle(lx, rx)
     end
-    return (lx + rx)/2
+    
+    return find_zero(x -> cdf(d, x) - p, (lx, rx), ITP(); xatol=tol)
 end
 
 function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::Real, rx::Real)
@@ -47,16 +32,14 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
-function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
+function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12, max_iters::Int=50)
     x = xs + (p - cdf(d, xs)) / pdf(d, xs)
     T = typeof(x)
+    F(x) = cdf(d, x) - p
+    f(x) = pdf(d, x)
     if 0 < p < 1
-        x0 = T(xs)
-        while abs(x-x0) > max(abs(x),abs(x0)) * tol
-            x0 = x
-            x = x0 + (p - cdf(d, x0)) / pdf(d, x0)
-        end
-        return x
+        tol = max(tol, cbrt(eps(float(T)))^2)
+        return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters)
     elseif p == 0
         return T(minimum(d))
     elseif p == 1
@@ -66,16 +49,14 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
     end
 end
 
-function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12)
+function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12, max_iters::Int=50)
     x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
     T = typeof(x)
+    F(x) = ccdf(d, x) - p
+    f(x) = -pdf(d, x)
     if 0 < p < 1
-        x0 = T(xs)
-        while abs(x-x0) > max(abs(x),abs(x0)) * tol
-            x0 = x
-            x = x0 + (ccdf(d, x0)-p) / pdf(d, x0)
-        end
-        return x
+        tol = max(tol, cbrt(eps(float(T)))^2)
+        return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters)
     elseif p == 1
         return T(minimum(d))
     elseif p == 0
@@ -85,24 +66,17 @@ function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real
     end
 end
 
-function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
+function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12, max_iters::Int=50)
     T = typeof(lp - logpdf(d,xs))
+    F(x) = logcdf(d, x) - lp
+    f(x) = exp(logpdf(d, x) - logcdf(d, x))
     if -Inf < lp < 0
         x0 = T(xs)
-        if lp < logcdf(d,x0)
-            x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
-                x0 = x
-                x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0)))
-            end
-        else
-            x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0))*tol
-                x0 = x
-                x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
-            end
-        end
-        return x
+        x = lp < logcdf(d,x0) ?
+            x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0))) :
+            x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0)))
+        tol = max(tol, cbrt(eps(float(T)))^2)
+        return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters)
     elseif lp == -Inf
         return T(minimum(d))
     elseif lp == 0
@@ -112,24 +86,17 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
     end
 end
 
-function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
+function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12, max_iters::Int=50)
     T = typeof(lp - logpdf(d,xs))
+    F(x) = logccdf(d, x) - lp
+    f(x) = -exp(logpdf(d, x) - logccdf(d, x))
     if -Inf < lp < 0
         x0 = T(xs)
-        if lp < logccdf(d,x0)
-            x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
-                x0 = x
-                x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0)))
-            end
-        else
-            x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
-            while abs(x-x0) >= max(abs(x),abs(x0)) * tol
-                x0 = x
-                x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
-            end
-        end
-        return x
+        x = lp < logccdf(d,x0) ?
+            x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0))) :
+            x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0)))
+        tol = max(tol, cbrt(eps(float(T)))^2)
+        return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters)
     elseif lp == -Inf
         return T(maximum(d))
     elseif lp == 0
diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl
index 3d319db25..a06dad7e8 100644
--- a/src/univariate/continuous/chernoff.jl
+++ b/src/univariate/continuous/chernoff.jl
@@ -189,14 +189,7 @@ function quantile(d::Chernoff, tau::Real)
     # one good approximation of the quantiles can be computed using Normal(0.0, stdapprox) with stdapprox = 0.52
     stdapprox = 0.52
     dnorm = Normal(0.0, 1.0)
-    if tau < 0.001
-        return -newton(x -> tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, 1.0 - tau)*stdapprox)
-
-    end
-    if tau > 0.999
-        return newton(x -> 1.0 - tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox)
-    end
-    return newton(x -> ChernoffComputations._cdf(x) - tau, ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox)   # should consider replacing x-> construct for speed
+    return quantile_newton(d, tau, quantile(dnorm, tau)*stdapprox)
 end
 
 minimum(d::Chernoff) = -Inf
diff --git a/test/mixture.jl b/test/mixture.jl
index d2b2742eb..1bec69472 100644
--- a/test/mixture.jl
+++ b/test/mixture.jl
@@ -290,4 +290,37 @@ end
             end
         end
     end
+
+    # issues #1611, #1807 and #1869
+    @testset "quantile convergence regressions" begin
+        cases = (
+            (
+                "nearly equal component quantiles (#1869)",
+                MixtureModel([Normal(0, 1), Normal(eps(), 1)], [0.999, 0.001]),
+                0.001,
+            ),
+            (
+                "adjacent floating point bracket (#1611)",
+                MixtureModel(
+                    [
+                        Uniform{Float64}(-0.0001, 0.0001),
+                        LogNormal{Float64}(11.174347445936371, 1.6086247197750911),
+                    ],
+                    [0.8832, 0.1168],
+                ),
+                0.99,
+            ),
+            (
+                "high quantile of exponential mixture (#1807)",
+                MixtureModel(Exponential, [10_000, 1_000_000], [0.5, 0.5]),
+                0.999,
+            ),
+        )
+
+        @testset "$name" for (name, d, p) in cases
+            x = quantile(d, p)
+            @test isfinite(x)
+            @test isapprox(cdf(d, x), p; atol=1e-10)
+        end
+    end
 end
diff --git a/test/quantile_newton.jl b/test/quantile_newton.jl
index da9291b0b..3b711fb26 100644
--- a/test/quantile_newton.jl
+++ b/test/quantile_newton.jl
@@ -5,3 +5,17 @@ using Test
 d = Normal()
 @test Distributions.quantile_newton(d, 0.5) == quantile(d, 0.5)
 @test Distributions.cquantile_newton(d, 0.5) == cquantile(d, 0.5)
+
+# issues #1571 and #2061
+@testset "InverseGaussian quantile convergence" begin
+    cases = (
+        (InverseGaussian(1.0, 0.25), 0.999996),
+        (InverseGaussian(2.8853900817779268), 0.9999996485182184),
+    )
+
+    @testset "p=$p" for (d, p) in cases
+        x = quantile(d, p)
+        @test isfinite(x)
+        @test isapprox(cdf(d, x), p; atol=1e-10)
+    end
+end
diff --git a/test/univariate/continuous/chernoff.jl b/test/univariate/continuous/chernoff.jl
index a5f28d02b..763a2ebe1 100644
--- a/test/univariate/continuous/chernoff.jl
+++ b/test/univariate/continuous/chernoff.jl
@@ -155,4 +155,15 @@
     for i=1:size(pdftest, 1)
         @test isapprox(pdf(d, pdftest[i, 1]), pdftest[i, 2] ; atol = 1e-6)
     end
+
+    # issue #1999
+    @testset "truncated quantile for non-precomputed probability (#1999)" begin
+        td = truncated(Chernoff(); lower=0.1)
+        p = 0.75
+        x = quantile(td, p)
+
+        @test isfinite(x)
+        @test x >= 0.1
+        @test isapprox(cdf(td, x), p; atol=1e-10)
+    end
 end