From ead5bf673ab744b9bf0a7cf5a291d17b4a2e30b0 Mon Sep 17 00:00:00 2001 From: cedri Date: Tue, 9 Jun 2026 16:06:13 -0400 Subject: [PATCH 1/9] Fixed the issue by adding a fallback to bissection --- src/quantilealgs.jl | 70 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 60 insertions(+), 10 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 8aa9f1b89a..22d825b3e6 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -47,12 +47,21 @@ 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, maxiter::Int=10_000) x = xs + (p - cdf(d, xs)) / pdf(d, xs) T = typeof(x) if 0 < p < 1 x0 = T(xs) - while abs(x-x0) > max(abs(x),abs(x0)) * tol + x2 = x0 + have_x2 = false + for _ in 1:maxiter + abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x + if have_x2 && isfinite(x) && isfinite(x2) && + abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + return quantile_bisect(d, p, minmax(x, x0)...) + end + x2 = x0 + have_x2 = true x0 = x x = x0 + (p - cdf(d, x0)) / pdf(d, x0) end @@ -66,12 +75,21 @@ 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, maxiter::Int=10_000) x = xs + (ccdf(d, xs)-p) / pdf(d, xs) T = typeof(x) if 0 < p < 1 x0 = T(xs) - while abs(x-x0) > max(abs(x),abs(x0)) * tol + x2 = x0 + have_x2 = false + for _ in 1:maxiter + abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x + if have_x2 && isfinite(x) && isfinite(x2) && + abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + return quantile_bisect(d, 1 - p, minmax(x, x0)...) + end + x2 = x0 + have_x2 = true x0 = x x = x0 + (ccdf(d, x0)-p) / pdf(d, x0) end @@ -85,19 +103,35 @@ 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, maxiter::Int=10_000) T = typeof(lp - logpdf(d,xs)) if -Inf < lp < 0 x0 = T(xs) + x2 = x0 + have_x2 = false 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 + for _ in 1:maxiter + abs(x-x0) < max(abs(x),abs(x0)) * tol && return x + if have_x2 && isfinite(x) && isfinite(x2) && + abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + return quantile_bisect(d, exp(lp), minmax(x, x0)...) + end + x2 = x0 + have_x2 = true 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 + for _ in 1:maxiter + abs(x-x0) < max(abs(x),abs(x0))*tol && return x + if have_x2 && isfinite(x) && isfinite(x2) && + abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + return quantile_bisect(d, exp(lp), minmax(x, x0)...) + end + x2 = x0 + have_x2 = true x0 = x x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0))) end @@ -112,19 +146,35 @@ 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, maxiter::Int=10_000) T = typeof(lp - logpdf(d,xs)) if -Inf < lp < 0 x0 = T(xs) + x2 = x0 + have_x2 = false 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 + for _ in 1:maxiter + abs(x-x0) < max(abs(x),abs(x0)) * tol && return x + if have_x2 && isfinite(x) && isfinite(x2) && + abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) + end + x2 = x0 + have_x2 = true 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 + for _ in 1:maxiter + abs(x-x0) < max(abs(x),abs(x0)) * tol && return x + if have_x2 && isfinite(x) && isfinite(x2) && + abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) + end + x2 = x0 + have_x2 = true x0 = x x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0))) end From 5f125dbd98178e510994f795665a9d3433b53d8d Mon Sep 17 00:00:00 2001 From: cedri Date: Mon, 15 Jun 2026 00:24:44 -0400 Subject: [PATCH 2/9] Improve quantile_newton stopping criterion - Add magnitude-aware tolerance floor to prevent chasing floating-point granularity - Reduce max_iter to 50 (proper stopping rule should rarely need it) - Add clear warnings for non-convergence with fallback to bisection - Detect 2-cycles explicitly as sign of ill-conditioning --- src/quantilealgs.jl | 68 +++++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 30 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 22d825b3e6..c12c69704d 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -47,17 +47,18 @@ 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, maxiter::Int=10_000) +function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=150) x = xs + (p - cdf(d, xs)) / pdf(d, xs) T = typeof(x) if 0 < p < 1 x0 = T(xs) x2 = x0 have_x2 = false - for _ in 1:maxiter + tol = max(tol, cbrt(eps(float(T)))^2) + for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && isfinite(x) && isfinite(x2) && - abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 + @warn "quantile_newton: 2-cycle detected, falling back to bisection" p=p maxlog=1 return quantile_bisect(d, p, minmax(x, x0)...) end x2 = x0 @@ -65,7 +66,8 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real= x0 = x x = x0 + (p - cdf(d, x0)) / pdf(d, x0) end - return x + @warn "quantile_newton: maximum iterations reached, falling back to bisection" p=p maxlog=1 + return quantile_bisect(d, p, minmax(x, x0)...) elseif p == 0 return T(minimum(d)) elseif p == 1 @@ -75,17 +77,18 @@ 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, maxiter::Int=10_000) +function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=150) x = xs + (ccdf(d, xs)-p) / pdf(d, xs) T = typeof(x) if 0 < p < 1 x0 = T(xs) x2 = x0 have_x2 = false - for _ in 1:maxiter + tol = max(tol, cbrt(eps(float(T)))^2) + for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && isfinite(x) && isfinite(x2) && - abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 + @warn "cquantile_newton: 2-cycle detected, falling back to bisection" p=p maxlog=1 return quantile_bisect(d, 1 - p, minmax(x, x0)...) end x2 = x0 @@ -93,7 +96,8 @@ function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real x0 = x x = x0 + (ccdf(d, x0)-p) / pdf(d, x0) end - return x + @warn "cquantile_newton: maximum iterations reached, falling back to bisection" p=p maxlog=1 + return quantile_bisect(d, 1 - p, minmax(x, x0)...) elseif p == 1 return T(minimum(d)) elseif p == 0 @@ -103,18 +107,19 @@ 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, maxiter::Int=10_000) +function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=150) T = typeof(lp - logpdf(d,xs)) if -Inf < lp < 0 x0 = T(xs) x2 = x0 have_x2 = false + tol = max(tol, cbrt(eps(float(T)))^2) if lp < logcdf(d,x0) x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0))) - for _ in 1:maxiter - abs(x-x0) < max(abs(x),abs(x0)) * tol && return x - if have_x2 && isfinite(x) && isfinite(x2) && - abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + for _ in 1:max_iter + abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x + if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 + @warn "invlogcdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 return quantile_bisect(d, exp(lp), minmax(x, x0)...) end x2 = x0 @@ -124,10 +129,10 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea end else x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0))) - for _ in 1:maxiter - abs(x-x0) < max(abs(x),abs(x0))*tol && return x - if have_x2 && isfinite(x) && isfinite(x2) && - abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + for _ in 1:max_iter + abs(x-x0) <= max(abs(x),abs(x0))*tol && return x + if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 + @warn "invlogcdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 return quantile_bisect(d, exp(lp), minmax(x, x0)...) end x2 = x0 @@ -136,7 +141,8 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0))) end end - return x + @warn "invlogcdf_newton: maximum iterations reached, falling back to bisection" lp=lp maxlog=1 + return quantile_bisect(d, exp(lp), minmax(x, x0)...) elseif lp == -Inf return T(minimum(d)) elseif lp == 0 @@ -146,18 +152,19 @@ 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, maxiter::Int=10_000) +function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=150) T = typeof(lp - logpdf(d,xs)) if -Inf < lp < 0 x0 = T(xs) x2 = x0 have_x2 = false + tol = max(tol, cbrt(eps(float(T)))^2) if lp < logccdf(d,x0) x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0))) - for _ in 1:maxiter - abs(x-x0) < max(abs(x),abs(x0)) * tol && return x - if have_x2 && isfinite(x) && isfinite(x2) && - abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + for _ in 1:max_iter + abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x + if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 + @warn "invlogccdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) end x2 = x0 @@ -167,10 +174,10 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re end else x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0))) - for _ in 1:maxiter - abs(x-x0) < max(abs(x),abs(x0)) * tol && return x - if have_x2 && isfinite(x) && isfinite(x2) && - abs(x - x2) <= max(abs(x), abs(x2)) * cbrt(eps(float(T)))^2 + for _ in 1:max_iter + abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x + if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 + @warn "invlogccdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) end x2 = x0 @@ -179,7 +186,8 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0))) end end - return x + @warn "invlogccdf_newton: maximum iterations reached, falling back to bisection" lp=lp maxlog=1 + return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) elseif lp == -Inf return T(maximum(d)) elseif lp == 0 From e09c958a9db2dbe71f3591d7129ba8184feeca4d Mon Sep 17 00:00:00 2001 From: cedri Date: Mon, 15 Jun 2026 08:07:23 -0400 Subject: [PATCH 3/9] Now uses ITP from Roots.jl for Newton quantile fallback --- Project.toml | 2 ++ src/Distributions.jl | 1 + src/quantilealgs.jl | 48 ++++++++++++++++++++++---------------------- 3 files changed, 27 insertions(+), 24 deletions(-) diff --git a/Project.toml b/Project.toml index f674a61855..18bde8fb6c 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" @@ -49,6 +50,7 @@ PDMats = "0.11.35" Printf = "<0.0.1, 1" QuadGK = "2" Random = "<0.0.1, 1" +Roots = "2" SparseArrays = "<0.0.1, 1" SpecialFunctions = "1.2, 2" StableRNGs = "1" diff --git a/src/Distributions.jl b/src/Distributions.jl index ed2a2cdde7..55121af773 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 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 c12c69704d..4430a623e2 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -47,7 +47,7 @@ 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, max_iter::Int=150) +function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=50) x = xs + (p - cdf(d, xs)) / pdf(d, xs) T = typeof(x) if 0 < p < 1 @@ -58,16 +58,16 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real= for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "quantile_newton: 2-cycle detected, falling back to bisection" p=p maxlog=1 - return quantile_bisect(d, p, minmax(x, x0)...) + @warn "quantile_newton: 2-cycle detected, falling back to ITP" p=p maxlog=1 + return find_zero(x -> cdf(d, x) - p, minmax(x, x0), ITP()) end x2 = x0 have_x2 = true x0 = x x = x0 + (p - cdf(d, x0)) / pdf(d, x0) end - @warn "quantile_newton: maximum iterations reached, falling back to bisection" p=p maxlog=1 - return quantile_bisect(d, p, minmax(x, x0)...) + @warn "quantile_newton: maximum iterations reached, falling back to ITP" p=p maxlog=1 + return find_zero(x -> cdf(d, x) - p, minmax(x, x0), ITP()) elseif p == 0 return T(minimum(d)) elseif p == 1 @@ -77,7 +77,7 @@ 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, max_iter::Int=150) +function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=50) x = xs + (ccdf(d, xs)-p) / pdf(d, xs) T = typeof(x) if 0 < p < 1 @@ -88,16 +88,16 @@ function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "cquantile_newton: 2-cycle detected, falling back to bisection" p=p maxlog=1 - return quantile_bisect(d, 1 - p, minmax(x, x0)...) + @warn "cquantile_newton: 2-cycle detected, falling back to ITP" p=p maxlog=1 + return find_zero(x -> ccdf(d, x) - p, minmax(x, x0), ITP()) end x2 = x0 have_x2 = true x0 = x x = x0 + (ccdf(d, x0)-p) / pdf(d, x0) end - @warn "cquantile_newton: maximum iterations reached, falling back to bisection" p=p maxlog=1 - return quantile_bisect(d, 1 - p, minmax(x, x0)...) + @warn "cquantile_newton: maximum iterations reached, falling back to ITP" p=p maxlog=1 + return find_zero(x -> ccdf(d, x) - p, minmax(x, x0), ITP()) elseif p == 1 return T(minimum(d)) elseif p == 0 @@ -107,7 +107,7 @@ 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, max_iter::Int=150) +function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=50) T = typeof(lp - logpdf(d,xs)) if -Inf < lp < 0 x0 = T(xs) @@ -119,8 +119,8 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogcdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 - return quantile_bisect(d, exp(lp), minmax(x, x0)...) + @warn "invlogcdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 + return find_zero(x -> logcdf(d, x) - lp, minmax(x, x0), ITP()) end x2 = x0 have_x2 = true @@ -132,8 +132,8 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0))*tol && return x if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogcdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 - return quantile_bisect(d, exp(lp), minmax(x, x0)...) + @warn "invlogcdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 + return find_zero(x -> logcdf(d, x) - lp, minmax(x, x0), ITP()) end x2 = x0 have_x2 = true @@ -141,8 +141,8 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0))) end end - @warn "invlogcdf_newton: maximum iterations reached, falling back to bisection" lp=lp maxlog=1 - return quantile_bisect(d, exp(lp), minmax(x, x0)...) + @warn "invlogcdf_newton: maximum iterations reached, falling back to ITP" lp=lp maxlog=1 + return find_zero(x -> logcdf(d, x) - lp, minmax(x, x0), ITP()) elseif lp == -Inf return T(minimum(d)) elseif lp == 0 @@ -152,7 +152,7 @@ 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, max_iter::Int=150) +function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12, max_iter::Int=50) T = typeof(lp - logpdf(d,xs)) if -Inf < lp < 0 x0 = T(xs) @@ -164,8 +164,8 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogccdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 - return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) + @warn "invlogccdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 + return find_zero(x -> logccdf(d, x) - lp, minmax(x, x0), ITP()) end x2 = x0 have_x2 = true @@ -177,8 +177,8 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re for _ in 1:max_iter abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogccdf_newton: 2-cycle detected, falling back to bisection" lp=lp maxlog=1 - return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) + @warn "invlogccdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 + return find_zero(x -> logccdf(d, x) - lp, minmax(x, x0), ITP()) end x2 = x0 have_x2 = true @@ -186,8 +186,8 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0))) end end - @warn "invlogccdf_newton: maximum iterations reached, falling back to bisection" lp=lp maxlog=1 - return quantile_bisect(d, -expm1(lp), minmax(x, x0)...) + @warn "invlogccdf_newton: maximum iterations reached, falling back to ITP" lp=lp maxlog=1 + return find_zero(x -> logccdf(d, x) - lp, minmax(x, x0), ITP()) elseif lp == -Inf return T(maximum(d)) elseif lp == 0 From 7e265d488ef78b1fcb116a9a39beaad1ac77e515 Mon Sep 17 00:00:00 2001 From: cedri Date: Tue, 16 Jun 2026 02:17:59 -0400 Subject: [PATCH 4/9] Refactor quantile algorithms to use Roots.jl Delegates core root-finding logic in quantile_newton functions to Roots.find_zero using Newton(). Upgrades quantile_bisect to use the ITP() bracketing solver as suggested. Retains precision-based xrtol tolerance scaling. --- src/Distributions.jl | 2 +- src/quantilealgs.jl | 138 ++++++++----------------------------------- 2 files changed, 24 insertions(+), 116 deletions(-) diff --git a/src/Distributions.jl b/src/Distributions.jl index 55121af773..d00daf121c 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -4,7 +4,7 @@ using StatsBase, PDMats, StatsFuns, Statistics using StatsFuns: logtwo, invsqrt2, invsqrt2π import QuadGK: quadgk -using Roots: find_zero, ITP +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 4430a623e2..593b8b22c9 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -15,23 +15,7 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx # ≈ 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 - 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,27 +31,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, max_iter::Int=50) +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) - x2 = x0 - have_x2 = false tol = max(tol, cbrt(eps(float(T)))^2) - for _ in 1:max_iter - abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "quantile_newton: 2-cycle detected, falling back to ITP" p=p maxlog=1 - return find_zero(x -> cdf(d, x) - p, minmax(x, x0), ITP()) - end - x2 = x0 - have_x2 = true - x0 = x - x = x0 + (p - cdf(d, x0)) / pdf(d, x0) - end - @warn "quantile_newton: maximum iterations reached, falling back to ITP" p=p maxlog=1 - return find_zero(x -> cdf(d, x) - p, minmax(x, x0), ITP()) + return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters) elseif p == 0 return T(minimum(d)) elseif p == 1 @@ -77,27 +48,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, max_iter::Int=50) +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) - x2 = x0 - have_x2 = false tol = max(tol, cbrt(eps(float(T)))^2) - for _ in 1:max_iter - abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "cquantile_newton: 2-cycle detected, falling back to ITP" p=p maxlog=1 - return find_zero(x -> ccdf(d, x) - p, minmax(x, x0), ITP()) - end - x2 = x0 - have_x2 = true - x0 = x - x = x0 + (ccdf(d, x0)-p) / pdf(d, x0) - end - @warn "cquantile_newton: maximum iterations reached, falling back to ITP" p=p maxlog=1 - return find_zero(x -> ccdf(d, x) - p, minmax(x, x0), ITP()) + return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters) elseif p == 1 return T(minimum(d)) elseif p == 0 @@ -107,42 +65,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, max_iter::Int=50) +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) - x2 = x0 - have_x2 = false + 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) - if lp < logcdf(d,x0) - x = x0 - exp(lp - logpdf(d,x0) + logexpm1(max(logcdf(d,x0)-lp,0))) - for _ in 1:max_iter - abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogcdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 - return find_zero(x -> logcdf(d, x) - lp, minmax(x, x0), ITP()) - end - x2 = x0 - have_x2 = true - 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))) - for _ in 1:max_iter - abs(x-x0) <= max(abs(x),abs(x0))*tol && return x - if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogcdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 - return find_zero(x -> logcdf(d, x) - lp, minmax(x, x0), ITP()) - end - x2 = x0 - have_x2 = true - x0 = x - x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0))) - end - end - @warn "invlogcdf_newton: maximum iterations reached, falling back to ITP" lp=lp maxlog=1 - return find_zero(x -> logcdf(d, x) - lp, minmax(x, x0), ITP()) + return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters) elseif lp == -Inf return T(minimum(d)) elseif lp == 0 @@ -152,42 +85,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, max_iter::Int=50) +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) - x2 = x0 - have_x2 = false + 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) - if lp < logccdf(d,x0) - x = x0 + exp(lp - logpdf(d,x0) + logexpm1(max(logccdf(d,x0)-lp,0))) - for _ in 1:max_iter - abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogccdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 - return find_zero(x -> logccdf(d, x) - lp, minmax(x, x0), ITP()) - end - x2 = x0 - have_x2 = true - 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))) - for _ in 1:max_iter - abs(x-x0) <= max(abs(x),abs(x0)) * tol && return x - if have_x2 && abs(x-x2) <= max(abs(x),abs(x2)) * cbrt(eps(float(T)))^2 - @warn "invlogccdf_newton: 2-cycle detected, falling back to ITP" lp=lp maxlog=1 - return find_zero(x -> logccdf(d, x) - lp, minmax(x, x0), ITP()) - end - x2 = x0 - have_x2 = true - x0 = x - x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0))) - end - end - @warn "invlogccdf_newton: maximum iterations reached, falling back to ITP" lp=lp maxlog=1 - return find_zero(x -> logccdf(d, x) - lp, minmax(x, x0), ITP()) + return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters) elseif lp == -Inf return T(maximum(d)) elseif lp == 0 From 0cd0cf7e5fd3a6b38407e022fec726d40d019860 Mon Sep 17 00:00:00 2001 From: cedri Date: Thu, 18 Jun 2026 01:13:21 -0400 Subject: [PATCH 5/9] - `quantile_bisect` now uses isapprox to check whether lx and rx are equal - the `chernoff` distribution now uses `quantile_newton` for none precomputed quantiles - Added test cases for the following issues (#1571, #1611, #1807, #1869, #1999, #2061) --- src/quantilealgs.jl | 15 ++++++------ src/univariate/continuous/chernoff.jl | 9 +------ test/mixture.jl | 33 ++++++++++++++++++++++++++ test/quantile_newton.jl | 14 +++++++++++ test/univariate/continuous/chernoff.jl | 11 +++++++++ 5 files changed, 67 insertions(+), 15 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 593b8b22c9..cf83aee93f 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -3,18 +3,19 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx::T) where {T<:Real} rx < lx && throw(ArgumentError("empty bracketing interval [$lx, $rx]")) + # 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 + # 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 + if isapprox(lx, rx; atol=tol, rtol=tol) # 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) + return middle(lx, rx) 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 + return find_zero(x -> cdf(d, x) - p, (lx, rx), ITP(); xatol=tol) end diff --git a/src/univariate/continuous/chernoff.jl b/src/univariate/continuous/chernoff.jl index 3d319db252..a06dad7e81 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 d2b2742ebc..1bec694721 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 da9291b0b8..3b711fb268 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 a5f28d02bc..763a2ebe1e 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 From a1ac512b5b4461fe174e1eba2c7f0b5252fa2dae Mon Sep 17 00:00:00 2001 From: cedri Date: Fri, 19 Jun 2026 00:45:11 -0400 Subject: [PATCH 6/9] Fix quantile bisection convergence for collapsed brackets --- src/quantilealgs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index cf83aee93f..538b8a01a1 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -10,7 +10,7 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, 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 isapprox(lx, rx; atol=tol, rtol=tol) + if 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) From 1f4cc9babdb11b93c12b0b33efbcb7f5da9514d8 Mon Sep 17 00:00:00 2001 From: cedri Date: Fri, 19 Jun 2026 16:02:28 -0400 Subject: [PATCH 7/9] Handle equal quantile bisection endpoints --- src/quantilealgs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 538b8a01a1..8976cab075 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -10,7 +10,7 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, 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 <= tol * max(abs(lx), abs(rx)) + 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) From 0dacff04e5e4e1f01ef54acc543998ec8edb84a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Fri, 19 Jun 2026 23:40:45 +0200 Subject: [PATCH 8/9] Use Roots defaults with a bracketing fallback for quantile root-finding MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Building on #2066, simplify the quantile solvers to rely on Roots' default tolerances instead of hand-tuned options: - `quantile_bisect` delegates to `find_zero(g, (lx, rx), ITP())`. ITP's relative `xrtol = eps` converges even for brackets far from zero, where the previous absolute `cbrt(eps)^2` tolerance fell below the floating-point spacing and looped forever (#1611, #1807). - The four `_newton` functions delegate to `find_zero((F, f), x, Newton(), ITP())`. Roots switches to the ITP bracketing method as soon as Newton brackets the root, which resolves the oscillation/stalling seen for extreme quantiles (#1571, #1898, #2061). The `tol`/`max_iters` parameters and the explicit tolerance floors are dropped: Roots' Newton defaults combine an absolute (`xatol = eps`) and a relative (`xrtol = eps`) x-tolerance, so they converge both near `p ≈ 1` and for roots near zero without bespoke floors. `quantile_bisect` keeps only the exact `rx == lx` degenerate guard (#1501); the approximate "collapsed bracket" handling is dropped, as a non-bracketing interval is a separate concern from the solver choice. Tests use `isapprox`'s type-appropriate defaults rather than fixed tolerances, and are split one regression per testset. Adds coverage for #1898 and a `Float32` round-trip. Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 2 +- src/quantilealgs.jl | 40 ++++++++++----------- test/mixture.jl | 49 +++++++++++--------------- test/quantile_newton.jl | 42 +++++++++++++++++----- test/univariate/continuous/chernoff.jl | 2 +- 5 files changed, 75 insertions(+), 60 deletions(-) diff --git a/Project.toml b/Project.toml index bb91ff1f63..716a31946d 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ PDMats = "0.11.35" Printf = "<0.0.1, 1" QuadGK = "2" Random = "<0.0.1, 1" -Roots = "2" +Roots = "3" SparseArrays = "<0.0.1, 1" SparseConnectivityTracer = "1" SpecialFunctions = "1.2, 2" diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 8976cab075..6ab2092089 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -3,20 +3,19 @@ function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::T, rx::T) where {T<:Real} rx < lx && throw(ArgumentError("empty bracketing interval [$lx, $rx]")) - # 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 - # 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)) + 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, rx) + return middle(lx) end - - return find_zero(x -> cdf(d, x) - p, (lx, rx), ITP(); xatol=tol) + + # Rely on Roots' default tolerances for the bracketing solver. Unlike the previous + # hand-rolled bisection, ITP uses a relative `xrtol = eps`, so it converges even for + # brackets far from zero where an absolute tolerance is below the floating-point spacing + # (#1611). + return find_zero(x -> cdf(d, x) - p, (lx, rx), ITP()) end function quantile_bisect(d::ContinuousUnivariateDistribution, p::Real, lx::Real, rx::Real) @@ -32,14 +31,16 @@ 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, max_iters::Int=50) +function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d)) 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 - tol = max(tol, cbrt(eps(float(T)))^2) - return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters) + # Newton's method with an ITP bracketing fallback: Roots switches to ITP whenever the + # iteration brackets the root (sign change of `F`), which resolves the oscillation/stalling + # observed for extreme quantiles (#2061, #1898). Roots' default tolerances are used. + return find_zero((F, f), x, Newton(), ITP()) elseif p == 0 return T(minimum(d)) elseif p == 1 @@ -49,14 +50,13 @@ 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, max_iters::Int=50) +function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d)) 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 - tol = max(tol, cbrt(eps(float(T)))^2) - return find_zero((F, f), x, Newton(), ITP(); xrtol=tol, maxiters=max_iters) + return find_zero((F, f), x, Newton(), ITP()) elseif p == 1 return T(minimum(d)) elseif p == 0 @@ -66,7 +66,7 @@ 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, max_iters::Int=50) +function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d)) T = typeof(lp - logpdf(d,xs)) F(x) = logcdf(d, x) - lp f(x) = exp(logpdf(d, x) - logcdf(d, x)) @@ -75,8 +75,7 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea 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) + return find_zero((F, f), x, Newton(), ITP()) elseif lp == -Inf return T(minimum(d)) elseif lp == 0 @@ -86,7 +85,7 @@ 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, max_iters::Int=50) +function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d)) T = typeof(lp - logpdf(d,xs)) F(x) = logccdf(d, x) - lp f(x) = -exp(logpdf(d, x) - logccdf(d, x)) @@ -95,8 +94,7 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re 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) + return find_zero((F, f), x, Newton(), ITP()) elseif lp == -Inf return T(maximum(d)) elseif lp == 0 diff --git a/test/mixture.jl b/test/mixture.jl index 1bec694721..92493f1a03 100644 --- a/test/mixture.jl +++ b/test/mixture.jl @@ -291,36 +291,27 @@ 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, - ), + # issue #1611: bracket of adjacent floating point numbers far from zero + @testset "adjacent floating point bracket (#1611)" begin + d = MixtureModel( + [ + Uniform{Float64}(-0.0001, 0.0001), + LogNormal{Float64}(11.174347445936371, 1.6086247197750911), + ], + [0.8832, 0.1168], ) + p = 0.99 + x = quantile(d, p) + @test isfinite(x) + @test cdf(d, x) ≈ p + end - @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 + # issue #1807: high quantile of a mixture of exponentials with very different scales + @testset "high quantile of exponential mixture (#1807)" begin + d = MixtureModel(Exponential, [10_000, 1_000_000], [0.5, 0.5]) + p = 0.999 + x = quantile(d, p) + @test isfinite(x) + @test cdf(d, x) ≈ p end end diff --git a/test/quantile_newton.jl b/test/quantile_newton.jl index 3b711fb268..bb9bcc1687 100644 --- a/test/quantile_newton.jl +++ b/test/quantile_newton.jl @@ -6,16 +6,42 @@ 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), - ) +# issue #1571 +@testset "InverseGaussian quantile convergence (#1571)" begin + d = InverseGaussian(1.0, 0.25) + p = 0.999996 + x = quantile(d, p) + @test isfinite(x) + @test cdf(d, x) ≈ p +end + +# issue #2061 +@testset "InverseGaussian quantile convergence (#2061)" begin + d = InverseGaussian(2.8853900817779268) + p = 0.9999996485182184 + x = quantile(d, p) + @test isfinite(x) + @test cdf(d, x) ≈ p +end + +# issue #1898: large-σ InverseGaussian, `p` far in the tail. `2 * cdf(Normal(), 5) - 1` is the +# value `erf(5 / sqrt(2))` from the original report. +@testset "InverseGaussian quantile convergence (#1898)" begin + d = InverseGaussian(1.187997687788096, 60.467382225458564) + p = 2 * cdf(Normal(), 5) - 1 + x = quantile(d, p) + @test isfinite(x) + @test cdf(d, x) ≈ p +end - @testset "p=$p" for (d, p) in cases +# the Roots defaults are type-aware: `Float32` must converge (and preserve the element type) +# without the absolute tolerance floor that the hand-rolled loops needed to avoid an infinite loop. +@testset "Float32 quantile round-trip" begin + d = InverseGaussian(2.0f0, 3.0f0) + for p in (0.1f0, 0.5f0, 0.9f0, 0.999f0) x = quantile(d, p) + @test x isa Float32 @test isfinite(x) - @test isapprox(cdf(d, x), p; atol=1e-10) + @test cdf(d, x) ≈ p end end diff --git a/test/univariate/continuous/chernoff.jl b/test/univariate/continuous/chernoff.jl index 763a2ebe1e..c72a3c0493 100644 --- a/test/univariate/continuous/chernoff.jl +++ b/test/univariate/continuous/chernoff.jl @@ -164,6 +164,6 @@ @test isfinite(x) @test x >= 0.1 - @test isapprox(cdf(td, x), p; atol=1e-10) + @test cdf(td, x) ≈ p end end From 2561f2eade752ff3f0a1f0726d4da485eb5f946d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Thu, 25 Jun 2026 14:53:19 +0200 Subject: [PATCH 9/9] Bump version from 0.25.127 to 0.25.128 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 716a31946d..72e8e16959 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.127" +version = "0.25.128" [deps] AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"