Skip to content

Quantile: use Newton method from Roots.jl #1900

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 18 commits into
base: master
Choose a base branch
from
Open

Quantile: use Newton method from Roots.jl #1900

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
c443cbd
`quantile_newton()`: extract `df()`
LebedevRI Sep 12, 2024
4a40f79
`cquantile_newton()`: extract `df()`
LebedevRI Sep 12, 2024
c02e6f7
`quantile_newton()`: extract `newton()`
LebedevRI Sep 12, 2024
e0928ba
`cquantile_newton()`: refactor using `newton()`
LebedevRI Sep 12, 2024
8b1267e
`invlogcdf_newton()`/`invlogccdf_newton()`: use strict inequality
LebedevRI Sep 12, 2024
905075a
`invlogcdf_newton()`: extract `df(x)`
LebedevRI Sep 12, 2024
8d0ebe2
`invlogccdf_newton()`: extract `df(x)`
LebedevRI Sep 12, 2024
efe20ec
`invlogcdf_newton()`: refactor using `newton()`
LebedevRI Sep 12, 2024
db6bf0a
`invlogccdf_newton()`: refactor using `newton()`
LebedevRI Sep 12, 2024
b97bfca
`newton()`: more closely match text-book method, which performs subtr…
LebedevRI Sep 12, 2024
fd6017d
`newton()`: `f()` is more of a function to calculate the actual delta
LebedevRI Sep 12, 2024
f3533d1
`[c]quantile_newton()`: explicitly pass function and derivative into …
LebedevRI Sep 12, 2024
46f4a10
`invlog[c]cdf_newton()`: explicitly pass function and derivative into…
LebedevRI Sep 13, 2024
4d021a7
`[c]quantile_newton()`: sink negation into `f(x)`
LebedevRI Sep 12, 2024
8399224
`newton_impl()`: use `isapprox()`
LebedevRI Sep 13, 2024
4a66424
Be more specific that `tol` is really a `xrtol`
LebedevRI Sep 13, 2024
53c2171
`newton()`: defer to `Roots.jl` implementation
LebedevRI Sep 12, 2024
f095427
Merge branch 'master' into roots.jl
LebedevRI May 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -48,6 +49,7 @@ PDMats = "0.11.35"
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"
Expand Down
70 changes: 32 additions & 38 deletions src/quantilealgs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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}
Expand Down Expand Up @@ -47,16 +49,19 @@ 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)
x = xs + (p - cdf(d, xs)) / pdf(d, xs)
function newton((f,df), xs::T=mode(d), xrtol::Real=1e-12) where {T}
return find_zero((f,df), xs, Roots.Newton(), xatol=0, xrtol=xrtol, atol=0, rtol=eps(float(T)), maxiters=typemax(Int))
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Notably, it would be sufficient to pass a reasonable maxiters to fix the currently-present bug.

end

function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)
f(x) = cdf(d, x) - p
df(x) = pdf(d, x)
# FIXME: can this be expressed via `promote_type()`? Test coverage missing.
Δ(x) = f(x)/df(x)
x = xs - Δ(xs)
T = typeof(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
return newton((f, df), T(xs), xrtol)
elseif p == 0
return T(minimum(d))
elseif p == 1
Expand All @@ -66,16 +71,15 @@ 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)
x = xs + (ccdf(d, xs)-p) / pdf(d, xs)
function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)
f(x) = p - ccdf(d, x)
df(x) = pdf(d, x)
# FIXME: can this be expressed via `promote_type()`? Test coverage missing.
Δ(x) = f(x)/df(x)
x = xs - Δ(xs)
T = typeof(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
return newton((f, df), T(xs), xrtol)
elseif p == 1
return T(minimum(d))
elseif p == 0
Expand All @@ -85,22 +89,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), xrtol::Real=1e-12)
T = typeof(lp - logpdf(d,xs))
f_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
f_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
df(x::T) where {T} = T(1)
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
return newton((f_ver0,df), T(xs), xrtol)
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
return newton((f_ver1,df), T(xs), xrtol)
end
return x
elseif lp == -Inf
Expand All @@ -112,22 +111,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), xrtol::Real=1e-12)
T = typeof(lp - logpdf(d,xs))
f_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
f_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
df(x::T) where {T} = T(1)
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
return newton((f_ver0,df), T(xs), xrtol)
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
return newton((f_ver1,df), T(xs), xrtol)
end
return x
elseif lp == -Inf
Expand Down
Loading