Limit number of Newton iterations in beta_inc_inv #396
Conversation
smaller starting values.
Codecov Report
@@ Coverage Diff @@
## master #396 +/- ##
===========================================
- Coverage 91.99% 37.30% -54.69%
===========================================
Files 12 12
Lines 2809 2761 -48
===========================================
- Hits 2584 1030 -1554
- Misses 225 1731 +1506
Flags with carried forward coverage won't be shown. Click here to find out more.
Continue to review full report at Codecov.
|
| @@ -968,8 +971,8 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p) | |||
| sq = 1.0 | |||
| prev = 1.0 | |||
|
|
|||
| if x < 0.0001 | |||
| x = 0.0001 | |||
| if x < 1e-200 | |||
There was a problem hiding this comment.
Is there any specific reason for this choice?
There was a problem hiding this comment.
Implementation-wise, it seems the section here could be simplified to
x = clamp(x, 1e-200, 0.9999)| @@ -1026,7 +1029,7 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p) | |||
| end | |||
|
|
|||
| #check if current estimate is acceptable | |||
|
|
|||
| prev, acu, p_approx, x, tx | |||
There was a problem hiding this comment.
I assume this was for debugging?
| prev, acu, p_approx, x, tx |
| @@ -1001,7 +1004,7 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p) | |||
| acu = exp10(iex) | |||
|
|
|||
| #iterate | |||
| while true | |||
| for i in 1:maxiter | |||
There was a problem hiding this comment.
Maybe
| for i in 1:maxiter | |
| for _ in 1:maxiter |
since we don't use i?
| @@ -1039,6 +1042,9 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p) | |||
| x = tx | |||
| p_approx_prev = p_approx | |||
| end | |||
|
|
|||
| @debug "Newton iterations didn't converge in $maxiter iterations. The result might have reduced precision." | |||
There was a problem hiding this comment.
This probably deserves a higher logging level?
| @debug "Newton iterations didn't converge in $maxiter iterations. The result might have reduced precision." | |
| @warn "Newton iterations didn't converge in $maxiter iterations. The result might have reduced precision." |
There was a problem hiding this comment.
I deliberately chose @debug over @warn to avoid that we raise warnings in normal usage. These functions sometimes lose precision in corner cases and it's not something we generally warn about. Eventually, we should have a better and more granular way of providing this information. E.g. via error codes and that's why I added the @debug statement.
|
The error here is because one of the existing test cases need a lot of iterations to converge. The starting values are really bad. I've read the two papers behind this implementation and I'm not trying to figure out if the starting values can be improved in the case that currently fails. |
* 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>
|
The issue was fixed by #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>
* 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>
and allow much smaller starting values.
I checked with the original Fortran version and it also hangs on the test problem. However, I think for different reasons since their tolerance is much higher than ours and the problem here is that we can't get below the
2.2250738585071875e-308tolerance. I'm a little skeptical that such a small tolerance is generally workable but the Fortran source comments thatI put the limit at 30 which is completely arbitrary but with a good starting value, I believe it should be fine. I decided to raise a debug warning when the function hits the iteration max. The clamping away from zero and one in
SpecialFunctions.jl/src/beta_inc.jl
Lines 971 to 976 in ce58a13