From c8264761a3069ca869ed9dbb474591c8743c75a4 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 05:10:10 +0300 Subject: [PATCH 01/12] `quantile_newton()`: extract `df()` --- src/quantilealgs.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 8aa9f1b89a..0cf41064ff 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -48,13 +48,14 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) = # 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) + f(x) = (p - cdf(d, x)) / pdf(d, x) + x = xs + f(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) + x = x0 + f(x0) end return x elseif p == 0 From 672e46b9cec8a209d0fcaecd2f26f2442bf87024 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 05:37:44 +0300 Subject: [PATCH 02/12] `cquantile_newton()`: extract `df()` --- src/quantilealgs.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 0cf41064ff..a1235df1fc 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -68,13 +68,14 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real= 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) + f(x) = (ccdf(d, x)-p) / pdf(d, x) + x = xs + f(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) + x = x0 + f(x0) end return x elseif p == 1 From fc9077839af61492ec2c32bf6e99b712c69d6696 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 05:38:16 +0300 Subject: [PATCH 03/12] `quantile_newton()`: extract `newton()` --- src/quantilealgs.jl | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index a1235df1fc..79bdc4be7f 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -47,17 +47,24 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) = # Distribution, with Application to the Inverse Gaussian Distribution # http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf +function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T} + x = xs + f(xs) + @assert typeof(x) === T + x0 = T(xs) + while abs(x-x0) > max(abs(x),abs(x0)) * tol + x0 = x + x = x0 + f(x0) + end + return x +end + function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12) f(x) = (p - cdf(d, x)) / pdf(d, x) + # FIXME: can this be expressed via `promote_type()`? Test coverage missing. x = xs + f(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 + f(x0) - end - return x + return newton(f, T(xs), tol) elseif p == 0 return T(minimum(d)) elseif p == 1 From 330f8fa23b62a8c6015b3515867fa47ac25ebd33 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 05:38:33 +0300 Subject: [PATCH 04/12] `cquantile_newton()`: refactor using `newton()` --- src/quantilealgs.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 79bdc4be7f..60e5034c50 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -76,15 +76,11 @@ end function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12) f(x) = (ccdf(d, x)-p) / pdf(d, x) + # FIXME: can this be expressed via `promote_type()`? Test coverage missing. x = xs + f(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 + f(x0) - end - return x + return newton(f, T(xs), tol) elseif p == 1 return T(minimum(d)) elseif p == 0 From 1f9b7a3fb3a54a0c2f6f545ca5fa9aa64069d0ae Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 06:20:51 +0300 Subject: [PATCH 05/12] `invlogcdf_newton()`/`invlogccdf_newton()`: use strict inequality This is what the `quantile_newton()`/`cquantile_newton()` does, because otherwise they were able to end up in an endless loops, when the initial point and the mode are the same, see #666. I'm not sure this is needed here, but the next change is going to refactor them to use general `newton()`, which would make this change anyway, so unless we need the current behaviour, let's do this change explicitly. --- src/quantilealgs.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 60e5034c50..b1d1de3494 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -96,13 +96,13 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea 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 + 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 + 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 @@ -123,13 +123,13 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re 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 + 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 + 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 From 59a28310deee3d520a9cebe34cbf3dc3734a08cd Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 06:49:46 +0300 Subject: [PATCH 06/12] `invlogcdf_newton()`: extract `df(x)` --- src/quantilealgs.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index b1d1de3494..d46ba013e2 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -92,19 +92,21 @@ end function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12) T = typeof(lp - logpdf(d,xs)) + f_a(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0))) + f_b(x) = exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0))) 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))) + x = x0 + f_a(x0) 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))) + x = x0 + f_a(x0) end else - x = x0 + exp(lp - logpdf(d,x0) + log1mexp(min(logcdf(d,x0)-lp,0))) + x = x0 + f_b(x0) 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))) + x = x0 + f_b(x0) end end return x From 2590ac0e654f65a387094797c5f662cba4402815 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 06:56:31 +0300 Subject: [PATCH 07/12] `invlogccdf_newton()`: extract `df(x)` --- src/quantilealgs.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index d46ba013e2..e7995b0746 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -121,19 +121,21 @@ end function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12) T = typeof(lp - logpdf(d,xs)) + f_a(x) = exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0))) + f_b(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0))) 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))) + x = x0 + f_a(x0) 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))) + x = x0 + f_a(x0) end else - x = x0 - exp(lp - logpdf(d,x0) + log1mexp(min(logccdf(d,x0)-lp,0))) + x = x0 + f_b(x0) 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))) + x = x0 + f_b(x0) end end return x From 4e7f76c5f0b038f4ea259f582028acc4ebd54ed3 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 07:25:36 +0300 Subject: [PATCH 08/12] `invlogcdf_newton()`: refactor using `newton()` --- src/quantilealgs.jl | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index e7995b0746..59c294f4f7 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -97,17 +97,9 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea if -Inf < lp < 0 x0 = T(xs) if lp < logcdf(d,x0) - x = x0 + f_a(x0) - while abs(x-x0) > max(abs(x),abs(x0)) * tol - x0 = x - x = x0 + f_a(x0) - end + return newton(f_a, T(xs), tol) else - x = x0 + f_b(x0) - while abs(x-x0) > max(abs(x),abs(x0))*tol - x0 = x - x = x0 + f_b(x0) - end + return newton(f_b, T(xs), tol) end return x elseif lp == -Inf From a319a7adc4f05fec2b14c5138db6ead0d16169b7 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 07:25:43 +0300 Subject: [PATCH 09/12] `invlogccdf_newton()`: refactor using `newton()` --- src/quantilealgs.jl | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 59c294f4f7..b987b39aec 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -118,17 +118,9 @@ function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Re if -Inf < lp < 0 x0 = T(xs) if lp < logccdf(d,x0) - x = x0 + f_a(x0) - while abs(x-x0) > max(abs(x),abs(x0)) * tol - x0 = x - x = x0 + f_a(x0) - end + return newton(f_a, T(xs), tol) else - x = x0 + f_b(x0) - while abs(x-x0) > max(abs(x),abs(x0)) * tol - x0 = x - x = x0 + f_b(x0) - end + return newton(f_b, T(xs), tol) end return x elseif lp == -Inf From 5e621bd1f2e7e2a7f600a6e247fd35306da40ea4 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 07:47:55 +0300 Subject: [PATCH 10/12] `newton()`: extract `converged()` --- src/quantilealgs.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index b987b39aec..a1a20e7470 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -48,10 +48,11 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) = # http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T} + converged(a,b) = abs(a-b) <= max(abs(a),abs(b)) * tol x = xs + f(xs) @assert typeof(x) === T x0 = T(xs) - while abs(x-x0) > max(abs(x),abs(x0)) * tol + while !converged(x, x0) x0 = x x = x0 + f(x0) end From 33c81f252ced79ab5eeb0af21828f38ffab84995 Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 08:00:48 +0300 Subject: [PATCH 11/12] `newton()`: make logic to keep previous "roots" more straight-forward --- src/quantilealgs.jl | 36 +++++++++++++++++++++++++++++------- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index a1a20e7470..1b4c83659c 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -41,6 +41,27 @@ end quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) = quantile_bisect(d, p, minimum(d), maximum(d)) +mutable struct ValueAndPredecessor{T} + storage::Tuple{T,T} + + function ValueAndPredecessor{T}(storage::Tuple{T,T}) where {T} + return new(storage) + end + + function ValueAndPredecessor{T}() where {T} + storage = Tuple{T,T}((T(NaN),T(NaN))) + return new(storage) + end + + function Base.push!(huh::ValueAndPredecessor{T}, e::T) where {T} + huh.storage = Tuple{T,T}((huh.storage[2], e)) + end + + function Base.getindex(ValueAndPredecessor::ValueAndPredecessor{T}, delay::Int)::T where {T} + return ValueAndPredecessor.storage[length(ValueAndPredecessor.storage) + delay] + end +end + # if starting at mode, Newton is convergent for any unimodal continuous distribution, see: # Göknur Giner, Gordon K. Smyth (2014) # A Monotonically Convergent Newton Iteration for the Quantiles of any Unimodal @@ -49,14 +70,15 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) = function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T} converged(a,b) = abs(a-b) <= max(abs(a),abs(b)) * tol - x = xs + f(xs) - @assert typeof(x) === T - x0 = T(xs) - while !converged(x, x0) - x0 = x - x = x0 + f(x0) + x = ValueAndPredecessor{T}() + push!(x, xs) + push!(x, x[0] + f(x[0])) + while !converged(x[0], x[-1]) + df = f(x[0]) + r = x[0] + df + push!(x, r) end - return x + return x[0] end function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), tol::Real=1e-12) From b0bc7372683677b77cd319689102a05066e46b1a Mon Sep 17 00:00:00 2001 From: Roman Lebedev Date: Thu, 12 Sep 2024 08:07:41 +0300 Subject: [PATCH 12/12] `newton()`: introduce oscillation dampening ("Newton-Roman") --- src/quantilealgs.jl | 8 ++++++ test/univariate/continuous/inversegaussian.jl | 25 ++++++++++++++++++- 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl index 1b4c83659c..69fa03a46d 100644 --- a/src/quantilealgs.jl +++ b/src/quantilealgs.jl @@ -76,6 +76,14 @@ function newton(f, xs::T=mode(d), tol::Real=1e-12) where {T} while !converged(x[0], x[-1]) df = f(x[0]) r = x[0] + df + # Vanilla Newton algorithm is known to be prone to oscillation, + # where we, e.g., go from 24.0 to 42.0, back to 24.0 and so forth. + # We can detect this situation by checking for convergence between + # the newly-computed "root" and the "root" we had two steps before. + if converged(r, x[-1]) + # To recover from oscillation, let's use just half of the delta. + r = r - (df / 2) + end push!(x, r) end return x[0] diff --git a/test/univariate/continuous/inversegaussian.jl b/test/univariate/continuous/inversegaussian.jl index 6990cae710..d22904ff0a 100644 --- a/test/univariate/continuous/inversegaussian.jl +++ b/test/univariate/continuous/inversegaussian.jl @@ -15,4 +15,27 @@ @test (p + q) ≈ 1 end end -end \ No newline at end of file +end + +@testset "InverseGaussian quantile" begin + p(num_σ) = erf(num_σ/sqrt(2)) + + begin + dist = InverseGaussian{Float64}(1.187997687788096, 60.467382225458564) + @test quantile(dist, p(4)) ≈ 1.9990956368573651 + @test quantile(dist, p(5)) ≈ 2.295607340999747 + @test quantile(dist, p(6)) ≈ 2.6249349452113484 + end + + @test quantile(InverseGaussian{Float64}(17.84806245738152, 163.707062977564), 0.9999981908772995) ≈ 69.37000274656731 + + begin + dist = InverseGaussian(1.0, 0.25) + @test quantile(dist, 0.99) ≈ 9.90306205018232 + @test quantile(dist, 0.999) ≈ 21.253279722084798 + @test quantile(dist, 0.9999) ≈ 34.73673452136752 + @test quantile(dist, 0.99999) ≈ 49.446586395457985 + @test quantile(dist, 0.999996) ≈ 55.53114044452607 + @test quantile(dist, 0.999999) ≈ 64.92521558088777 + end +end