Add cdf/ccdf/quantile to SkewNormal via Owen's T#2069
Open
devmotion wants to merge 12 commits into
Open
Conversation
- Add magnitude-aware tolerance floor to prevent chasing floating-point granularity - Reduce max_iter to 50 (proper stopping rule should rarely need it) - Add clear warnings for non-convergence with fallback to bisection - Detect 2-cycles explicitly as sign of ill-conditioning
Delegates core root-finding logic in quantile_newton functions to Roots.find_zero using Newton(). Upgrades quantile_bisect to use the ITP() bracketing solver as suggested. Retains precision-based xrtol tolerance scaling.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #2069 +/- ##
==========================================
+ Coverage 86.74% 86.89% +0.14%
==========================================
Files 149 149
Lines 8882 8898 +16
==========================================
+ Hits 7705 7732 +27
+ Misses 1177 1166 -11 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Building on #2066, simplify the quantile solvers to rely on Roots' default tolerances instead of hand-tuned options: - `quantile_bisect` delegates to `find_zero(g, (lx, rx), ITP())`. ITP's relative `xrtol = eps` converges even for brackets far from zero, where the previous absolute `cbrt(eps)^2` tolerance fell below the floating-point spacing and looped forever (#1611, #1807). - The four `_newton` functions delegate to `find_zero((F, f), x, Newton(), ITP())`. Roots switches to the ITP bracketing method as soon as Newton brackets the root, which resolves the oscillation/stalling seen for extreme quantiles (#1571, #1898, #2061). The `tol`/`max_iters` parameters and the explicit tolerance floors are dropped: Roots' Newton defaults combine an absolute (`xatol = eps`) and a relative (`xrtol = eps`) x-tolerance, so they converge both near `p ≈ 1` and for roots near zero without bespoke floors. `quantile_bisect` keeps only the exact `rx == lx` degenerate guard (#1501); the approximate "collapsed bracket" handling is dropped, as a non-bracketing interval is a separate concern from the solver choice. Tests use `isapprox`'s type-appropriate defaults rather than fixed tolerances, and are split one regression per testset. Adds coverage for #1898 and a `Float32` round-trip. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Implement the previously-missing (c)cdf for SkewNormal, now that StatsFuns
2.2 provides Owen's T function `owens_t`. The closed forms are
cdf(x) = Φ(z) - 2·T(z, α)
ccdf(x) = Φ̄(z) + 2·T(z, α), z = (x - ξ)/ω
with `logcdf`/`logccdf` using the generic `log(cdf)`/`log(ccdf)` fallbacks.
`pdf`/`logpdf` are refactored for consistency (destructure once, reuse `z`).
`quantile`/`cquantile` are added via Newton iteration from the mode, using
`quantile_newton`/`cquantile_newton`. This branch is stacked on #2070, whose
Roots-based solver handles the near-zero-root convergence that SkewNormal's
support (straddling 0) requires, so no change to `src/quantilealgs.jl` is
needed here. `invlogcdf`/`invlogccdf` are left to the generic
`quantile(d, exp(lp))` fallback: the log-space `invlogcdf_newton` overshoots
into the underflowed tail from the approximate mode `m_0` and returns NaN
there, whereas the generic fallback through `quantile_newton` stays robust.
Require StatsFuns 2.2 and bump the package version.
Testing: add a `sn`-based reference class (test/ref/continuous/skewnormal.R)
to the R reference framework, pin `sn` and its deps in renv.lock, and add the
generated entries to continuous_test.ref.json. The dedicated testset is
reduced to construction, support boundaries and the α=0 reduction to Normal.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
c027052 to
46b0d8d
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
SkewNormalwas missingcdf/ccdf(the source had#cdf requires Owen's T function.). StatsFuns 2.2 now providesowens_t, which unblocks the closed form:This PR adds
cdf/ccdf, refactorspdf/logpdffor consistency, and addsquantile/cquantile.logcdf/logccdf/invlogcdf/invlogccdfuse the generic fallbacks.Quantile root-finding (handled by #2070)
quantile/cquantileare Newton iterations from the mode viaquantile_newton/cquantile_newton. The original version of this PR also patchedsrc/quantilealgs.jlwith an absolute-tolerance floor, because the hand-rolled relative-only stopping rule never converges as the root approaches zero (whichSkewNormalhits). That fix is now subsumed by #2070, which replaces the hand-rolled loops withRoots.find_zeroand so converges for near-zero roots out of the box — this PR therefore drops its ownquantilealgs.jlchange.invlogcdf/invlogccdfare intentionally left to the genericquantile(d, exp(lp))fallback rather than the@quantile_newtonmacro: the log-spaceinvlogcdf_newtonovershoots into the underflowed tail from the approximate modem_0and returnsNaNfor extreme log-probabilities (e.g.p ≲ 1e-8), even with #2070's bracketing fallback, whereas routing throughquantile_newtonstays robust there.Tests
sn-based reference classtest/ref/continuous/skewnormal.R(moments fromsn:::sn.cumulants, an independent computation rather than a reimplementation of the Julia formulas), 4 entries tocontinuous_test.lst, the generated entries tocontinuous_test.ref.json, and pinnedsn+ deps inrenv.lock.Normal.Verification
snand Julia agree to ~1e-14.quantile/cquantileround-trips and theinvlogcdftail fallback verified against Simplify quantile root-finding to Roots defaults #2070's solver (incl. roots straddling 0).🤖 Generated with Claude Code