@@ -47,31 +47,31 @@ quantile_bisect(d::ContinuousUnivariateDistribution, p::Real) =
47
47
# Distribution, with Application to the Inverse Gaussian Distribution
48
48
# http://www.statsci.org/smyth/pubs/qinvgaussPreprint.pdf
49
49
50
- function newton_impl (Δ, xs:: T = mode (d), tol :: Real = 1e-12 ) where {T}
50
+ function newton_impl (Δ, xs:: T = mode (d), xrtol :: Real = 1e-12 ) where {T}
51
51
x = xs - Δ (xs)
52
52
@assert typeof (x) === T
53
53
x0 = T (xs)
54
- while ! isapprox (x, x0, atol= 0 , rtol= tol )
54
+ while ! isapprox (x, x0, atol= 0 , rtol= xrtol )
55
55
x0 = x
56
56
x = x0 - Δ (x0)
57
57
end
58
58
return x
59
59
end
60
60
61
- function newton ((f,df), xs:: T = mode (d), tol :: Real = 1e-12 ) where {T}
61
+ function newton ((f,df), xs:: T = mode (d), xrtol :: Real = 1e-12 ) where {T}
62
62
Δ (x) = f (x)/ df (x)
63
- return newton_impl (Δ, xs, tol )
63
+ return newton_impl (Δ, xs, xrtol )
64
64
end
65
65
66
- function quantile_newton (d:: ContinuousUnivariateDistribution , p:: Real , xs:: Real = mode (d), tol :: Real = 1e-12 )
66
+ function quantile_newton (d:: ContinuousUnivariateDistribution , p:: Real , xs:: Real = mode (d), xrtol :: Real = 1e-12 )
67
67
f (x) = cdf (d, x) - p
68
68
df (x) = pdf (d, x)
69
69
# FIXME : can this be expressed via `promote_type()`? Test coverage missing.
70
70
Δ (x) = f (x)/ df (x)
71
71
x = xs - Δ (xs)
72
72
T = typeof (x)
73
73
if 0 < p < 1
74
- return newton ((f, df), T (xs), tol )
74
+ return newton ((f, df), T (xs), xrtol )
75
75
elseif p == 0
76
76
return T (minimum (d))
77
77
elseif p == 1
@@ -81,15 +81,15 @@ function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=
81
81
end
82
82
end
83
83
84
- function cquantile_newton (d:: ContinuousUnivariateDistribution , p:: Real , xs:: Real = mode (d), tol :: Real = 1e-12 )
84
+ function cquantile_newton (d:: ContinuousUnivariateDistribution , p:: Real , xs:: Real = mode (d), xrtol :: Real = 1e-12 )
85
85
f (x) = p - ccdf (d, x)
86
86
df (x) = pdf (d, x)
87
87
# FIXME : can this be expressed via `promote_type()`? Test coverage missing.
88
88
Δ (x) = f (x)/ df (x)
89
89
x = xs - Δ (xs)
90
90
T = typeof (x)
91
91
if 0 < p < 1
92
- return newton ((f, df), T (xs), tol )
92
+ return newton ((f, df), T (xs), xrtol )
93
93
elseif p == 1
94
94
return T (minimum (d))
95
95
elseif p == 0
@@ -99,17 +99,17 @@ function cquantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real
99
99
end
100
100
end
101
101
102
- function invlogcdf_newton (d:: ContinuousUnivariateDistribution , lp:: Real , xs:: Real = mode (d), tol :: Real = 1e-12 )
102
+ function invlogcdf_newton (d:: ContinuousUnivariateDistribution , lp:: Real , xs:: Real = mode (d), xrtol :: Real = 1e-12 )
103
103
T = typeof (lp - logpdf (d,xs))
104
104
f_ver0 (x) = exp (lp - logpdf (d,x) + logexpm1 (max (logcdf (d,x)- lp,0 )))
105
105
f_ver1 (x) = - exp (lp - logpdf (d,x) + log1mexp (min (logcdf (d,x)- lp,0 )))
106
106
df (x:: T ) where {T} = T (1 )
107
107
if - Inf < lp < 0
108
108
x0 = T (xs)
109
109
if lp < logcdf (d,x0)
110
- return newton ((f_ver0,df), T (xs), tol )
110
+ return newton ((f_ver0,df), T (xs), xrtol )
111
111
else
112
- return newton ((f_ver1,df), T (xs), tol )
112
+ return newton ((f_ver1,df), T (xs), xrtol )
113
113
end
114
114
return x
115
115
elseif lp == - Inf
@@ -121,17 +121,17 @@ function invlogcdf_newton(d::ContinuousUnivariateDistribution, lp::Real, xs::Rea
121
121
end
122
122
end
123
123
124
- function invlogccdf_newton (d:: ContinuousUnivariateDistribution , lp:: Real , xs:: Real = mode (d), tol :: Real = 1e-12 )
124
+ function invlogccdf_newton (d:: ContinuousUnivariateDistribution , lp:: Real , xs:: Real = mode (d), xrtol :: Real = 1e-12 )
125
125
T = typeof (lp - logpdf (d,xs))
126
126
f_ver0 (x) = - exp (lp - logpdf (d,x) + logexpm1 (max (logccdf (d,x)- lp,0 )))
127
127
f_ver1 (x) = exp (lp - logpdf (d,x) + log1mexp (min (logccdf (d,x)- lp,0 )))
128
128
df (x:: T ) where {T} = T (1 )
129
129
if - Inf < lp < 0
130
130
x0 = T (xs)
131
131
if lp < logccdf (d,x0)
132
- return newton ((f_ver0,df), T (xs), tol )
132
+ return newton ((f_ver0,df), T (xs), xrtol )
133
133
else
134
- return newton ((f_ver1,df), T (xs), tol )
134
+ return newton ((f_ver1,df), T (xs), xrtol )
135
135
end
136
136
return x
137
137
elseif lp == - Inf
0 commit comments