From b0bc7372683677b77cd319689102a05066e46b1a Mon Sep 17 00:00:00 2001
From: Roman Lebedev <lebedev.ri@gmail.com>
Date: Thu, 12 Sep 2024 08:07:41 +0300
Subject: [PATCH] `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 1b4c83659..69fa03a46 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 6990cae71..d22904ff0 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