newton(): more closely match text-book method, which performs subtraction,

Author: LebedevRI

src/quantilealgs.jl

@@ -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}
@@ -48,20 +50,20 @@ 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}
-    x = xs + f(xs)
+    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)
+        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)
+    f(x) = -((p - cdf(d, x)) / pdf(d, x))
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
-    x = xs + f(xs)
+    x = xs - f(xs)
     T = typeof(x)
     if 0 < p < 1
         return newton(f, T(xs), tol)
@@ -75,9 +77,9 @@ 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)
-    f(x) = (ccdf(d, x)-p) / pdf(d, x)
+    f(x) = -((ccdf(d, x)-p) / pdf(d, x))
     # FIXME: can this be expressed via `promote_type()`? Test coverage missing.
-    x = xs + f(xs)
+    x = xs - f(xs)
     T = typeof(x)
     if 0 < p < 1
         return newton(f, T(xs), tol)
@@ -92,8 +94,8 @@ 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)))
+    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)
@@ -113,8 +115,8 @@ 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)))
+    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)