newton(): introduce oscillation dampening ("Newton-Roman")

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]

test/univariate/continuous/inversegaussian.jl

@@ -15,4 +15,27 @@
             @test (p + q) ≈ 1
         end
     end
-end
+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