Some fixes and simplications to the beta function(s). #316
Conversation
|
Hm. So the assertions were triggered in the incomplete beta tests. I'll look into it. |
There was a problem hiding this comment.
Looks mostly fine. I had some questions regarding the use of @assert and the changes in logabsbeta. In general, I wonder if it would be beneficial to use Float64 literals in all comparisons where the involved variables are guaranteed to be of type Float64 (apart from consistency).
| @@ -13,6 +13,7 @@ loggammadiv(a::Number, b::Number) = _loggammadiv(float(a), float(b)) | |||
|
|
|||
| _loggammadiv(a::T, b::T) where {T<:Base.IEEEFloat} = T(_loggammadiv(Float64(a), Float64(b))) # handle Float16, Float32 | |||
| function _loggammadiv(a::Float64, b::Float64) | |||
| @assert b >= 8 | |||
There was a problem hiding this comment.
I usually try to avoid @assert and use error instead, mainly due to the warning in the docstring of @assert that the assert statement might be disabled at various optimization levels and that it should be used only as a debugging tool. However, I guess you are more familiar with these Julia internals and I assume it does not matter as long as the behaviour does not depend critically on the assert statements.
There was a problem hiding this comment.
I do think it's only for debugging here. These function are only supposed to be used as part of other functions and it will actually be useful that they can be removed by optimization passes.
src/beta_inc.jl
Outdated
| @@ -348,6 +355,9 @@ External links: [DLMF](https://dlmf.nist.gov/8.17.22), [Wikipedia](https://en.wi | |||
| See also: [`beta_inc`](@ref) | |||
| """ | |||
| function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::Float64, w::Float64, epps::Float64) | |||
| @assert a >= 15 | |||
| @assert b <= 1 | |||
There was a problem hiding this comment.
Again, the docstring seems to indicate that also b < epps is required (a >= 15 and hence min(epps, epps * a) = epps). I assume that usually epps <= 1 and thus
| @assert b <= 1 | |
| @assert b < epps |
would be sufficient but probably safest would be
| @assert b <= 1 | |
| @assert b <= 1.0 && b < epps |
(b < min(1, epps) would be "incorrect" if epps > 1 and b <= min(1, epps) would be "incorrect" if epps < 1).
There was a problem hiding this comment.
This was the one that triggered the test failure so I had to check the Fortran source. First part of the comment here was simply wrong and I believe the assertion are correct.
| # asymptotic expansion for large `b` and positive arguments | ||
| # FIXME! We probably lose precision for negative arguments. However, `loggammadiv` | ||
| # is based on the NSWC code which doesn't bother with negative arguments | ||
| if a > 0 && b > 8 |
There was a problem hiding this comment.
Are these thresholds based on the NSWC code?
There was a problem hiding this comment.
Yes. All the beta functions in NSWC assume positive arguments and the cutoff used there is b>8.
of large magnitude differences was only applied to logabsbeta but not logbeta. They now share code so the special handling applies to both functions. The cutoff for special handling has also been lowered to avoid loss of precision for intermediate argument size differences. Finally, helper functions for the beta function now have assertions to ensure that they are only applied for the ranges they are written for.
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Codecov Report
@@ Coverage Diff @@
## master #316 +/- ##
==========================================
+ Coverage 89.28% 91.90% +2.61%
==========================================
Files 12 12
Lines 2633 2645 +12
==========================================
+ Hits 2351 2431 +80
+ Misses 282 214 -68
Flags with carried forward coverage won't be shown. Click here to find out more.
Continue to review full report at Codecov.
|
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
|
I believe I've addressed your comments but please let me know if I've missed anything. (I've added a few more tests to increase to coverage of some of the helper functions). |
direction is finite.
Float64 arguments.
|
Bump. I've just added a commit that changes some of the comparisons to be using float literals. |
| x = aa/ (aa+bb*exp(w^2)) | ||
| h = 2.0/(s + t) | ||
| w = pp_approx*sqrt(h + r)/h - (t - s)*(r + 5.0/6.0 - 2.0/(3.0*h)) | ||
| x = aa/(aa + bb*exp(w + w)) |
There was a problem hiding this comment.
So the w^2 in the original code was a typo? exp(w^2) ≠ exp(w + w)
There was a problem hiding this comment.
It is a typo. This is the line from the paper
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
The handling of large magnitude differences was only applied to logabsbeta but
not logbeta. They now share code so the special handling applies to
both functions. The cutoff for special handling has also been lowered
to avoid loss of precision for intermediate argument size differences.
Finally, helper functions for the beta function now have assertions
to ensure that they are only applied for the ranges they are written for.
@devmotion Would be great if you could review this one as well. The issue showed up when I tried to use native Julia implementations for the beta and binomial distributions in JuliaStats/StatsFuns.jl#113.