Skip to content

Commit 46f4a10

Browse files
committed
invlog[c]cdf_newton(): explicitly pass function and derivative into the `newton()
These have `df(x)=1`, at least that is what the `df(x)` is if we solve `Δ(x)=f(x)/f'(x)` as an equation for `f'(x)`.
1 parent f3533d1 commit 46f4a10

File tree

1 file changed

+10
-8
lines changed

1 file changed

+10
-8
lines changed

src/quantilealgs.jl

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -101,14 +101,15 @@ end
101101

102102
function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
103103
T = typeof(lp - logpdf(d,xs))
104-
Δ_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
105-
Δ_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
104+
f_ver0(x) = exp(lp - logpdf(d,x) + logexpm1(max(logcdf(d,x)-lp,0)))
105+
f_ver1(x) = -exp(lp - logpdf(d,x) + log1mexp(min(logcdf(d,x)-lp,0)))
106+
df(x::T) where {T} = T(1)
106107
if -Inf < lp < 0
107108
x0 = T(xs)
108109
if lp < logcdf(d,x0)
109-
return newton_impl(Δ_ver0, T(xs), tol)
110+
return newton((f_ver0,df), T(xs), tol)
110111
else
111-
return newton_impl(Δ_ver1, T(xs), tol)
112+
return newton((f_ver1,df), T(xs), tol)
112113
end
113114
return x
114115
elseif lp == -Inf
@@ -122,14 +123,15 @@ end
122123

123124
function invlogccdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Real=mode(d), tol::Real=1e-12)
124125
T = typeof(lp - logpdf(d,xs))
125-
Δ_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
126-
Δ_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
126+
f_ver0(x) = -exp(lp - logpdf(d,x) + logexpm1(max(logccdf(d,x)-lp,0)))
127+
f_ver1(x) = exp(lp - logpdf(d,x) + log1mexp(min(logccdf(d,x)-lp,0)))
128+
df(x::T) where {T} = T(1)
127129
if -Inf < lp < 0
128130
x0 = T(xs)
129131
if lp < logccdf(d,x0)
130-
return newton_impl(Δ_ver0, T(xs), tol)
132+
return newton((f_ver0,df), T(xs), tol)
131133
else
132-
return newton_impl(Δ_ver1, T(xs), tol)
134+
return newton((f_ver1,df), T(xs), tol)
133135
end
134136
return x
135137
elseif lp == -Inf

0 commit comments

Comments
 (0)