Fix hangs in beta_inc_inv (alternative/extension of #396)
#399
Conversation
| @@ -934,7 +934,6 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p) | |||
| r = sqrt(-2*log(p)) | |||
| p_approx = r - @horner(r, 2.30753e+00, 0.27061e+00) / @horner(r, 1.0, .99229e+00, .04481e+00) | |||
| fpu = floatmin(Float64) | |||
| sae = log10(fpu) | |||
There was a problem hiding this comment.
sae is only used in the computation of acu = exp10(max(..., sae)) below. I rewrote this as acu = max(exp10(...), fpu) (or rather 10 * fpu).
Codecov Report
@@ Coverage Diff @@
## master #399 +/- ##
==========================================
+ Coverage 91.99% 93.57% +1.58%
==========================================
Files 12 12
Lines 2809 2769 -40
==========================================
+ Hits 2584 2591 +7
+ Misses 225 178 -47
Flags with carried forward coverage won't be shown. Click here to find out more.
Continue to review full report at Codecov.
|
| g = 1.0 | ||
|
|
||
| tx = x - g*p_approx | ||
| while true | ||
| adj = g*p_approx | ||
| sq = adj^2 | ||
| adj = p_approx | ||
| tx = x - adj | ||
| while prev <= (sq = adj^2) || tx < 0.0 || tx > 1.0 | ||
| adj /= 3.0 | ||
| tx = x - adj | ||
| if (prev > sq && tx >= 0.0 && tx <= 1.0) | ||
| break | ||
| end | ||
| g /= 3.0 |
There was a problem hiding this comment.
I simplified this loop. Currently,
- the initial
adj = g*p_approxandtxis computed twice adjis recomputed in every step asg*p_approxwheregis divided by 3 in every iteration
In this PR, instead
- values in the initial iteration are only computed once
adjis scaled directly, without variablegand additional multiplication
| iex = max(-5.0/a^2 - 1.0/p^0.2 - 34.0, sae) | ||
| acu = exp10(iex) | ||
| iex = -5.0/a^2 - 1.0/p^0.2 - 34.0 | ||
| acu = max(exp10(iex), 10 * fpu) # 10 * fpu instead of fpu avoids hangs for small values |
There was a problem hiding this comment.
I think it's okay. Generally, I believe this calculation is pretty ad hoc anyway.
* Limit number of Newton iterations in beta_inc_inv and allow much smaller starting values. * Adjust tolerance instead of number of iterations Co-authored-by: Andreas Noack <andreas@noack.dk>
* Fix hangs in `beta_inc_inv` (alternative/extension of #396) (#399) * Limit number of Newton iterations in beta_inc_inv and allow much smaller starting values. * Adjust tolerance instead of number of iterations Co-authored-by: Andreas Noack <andreas@noack.dk> * Release 1.8.5 Co-authored-by: Andreas Noack <andreas@noack.dk>
| if x > .9999 | ||
| x = .9999 | ||
| end | ||
| x = clamp(x, fpu, prevfloat(1.0)) |
There was a problem hiding this comment.
This change, or rather the lower bound fpu, break examples such as JuliaStats/StatsFuns.jl#145. Unfortunately, the alternative 1e-200 proposed in #396 breaks that example as well.
I noticed that #396 is not merged yet, due to some minor comments and failing tests (I assume).
This PR builds on #396 and in addition to some minor changes (removal of line for debugging,
clampinstead of explicitifstatements) most importantly it uses a different approach for fixing the original issue: Lower-bounding the tolerance at10 * floatmin(Float64)instead offloatmin(Float64)seems to fix the pathological example.Of course, as the maximum number of iterations in #396 this is also an arbitrary choice. However, in contrast to the other PR tests still pass since harder examples that require a larger number of iterations don't fail.