newton(): make logic to keep previous "roots" more straight-forward

Author: LebedevRI

src/quantilealgs.jl

@@ -47,16 +47,44 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
 #   Distribution, with Application to the Inverse Gaussian Distribution
 #   http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
 
+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
+#   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}
     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)