Implement beta/gamma invlogcdf and invlogccdf in pure Julia#204
Implement beta/gamma invlogcdf and invlogccdf in pure Julia#204andreasnoack wants to merge 1 commit into
Conversation
andreasnoack
force-pushed
the
an/normath-betagamma
branch
from
c12b740 to
d7c7e63
Compare
Replace the Rmath-dependent invlogcdf/invlogccdf functions for beta and gamma distributions with Newton iteration in log-space. Uses asymptotic initial estimates for extreme tails and -expm1(lp) for accuracy near lp=0. Leverages existing betainvcdf/gammainvcdf for moderate cases.
andreasnoack
force-pushed
the
an/normath-betagamma
branch
from
d7c7e63 to
98c7729
Compare
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #204 +/- ##
==========================================
+ Coverage 72.39% 73.76% +1.37%
==========================================
Files 22 22
Lines 1007 1117 +110
==========================================
+ Hits 729 824 +95
- Misses 278 293 +15 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
| _lp = Float64(lp) | ||
|
|
||
| if isnan(_lp) || _lp > 0 | ||
| return convert(T, NaN) | ||
| elseif _lp == 0 | ||
| return one(T) | ||
| elseif isinf(_lp) # -Inf | ||
| return zero(T) | ||
| end | ||
|
|
||
| # Handle degenerate cases | ||
| _α = Float64(α); _β = Float64(β) |
There was a problem hiding this comment.
These Float64 calls will error for Duals and (arguably even worse) lead to silently incorrect results for higher-precision numbers. Ideally, we should just avoid the Float64 conversions in this function and instead ensure that the functions called here are written generically enough to support all common use cases. Then the responsibility for potentially handling other or higher-precision Reals would be off-loaded to those.
There was a problem hiding this comment.
I should probably convert the PRs to drafts. The current ones are mostly how Claude Code prepared them and I need to go through the code in detail.
I agree that a function like this should probably avoid conversion to Float64, but most of the functions containing the actual algorithms will necessarily be restricted to Float64 because the number of terms in the various approximations has been determined based on the precision.
3 participants
Summary
betainvlogcdf,betainvlogccdf,gammainvlogcdf,gammainvlogccdfwith pure Julia implementationsbetainvcdf/gammainvcdfvia-expm1(lp)for accuracy nearlp=0exp(lp)underflows): Newton iteration in log-space usingbetalogcdf/gammalogcdfandbetalogpdf/gammalogpdf, with asymptotic initial estimatesTest plan
rmathcomp_tests)lp = -5000match Rmath where representable🤖 Generated with Claude Code