Skip to content

Add cdf/ccdf/quantile to SkewNormal via Owen's T#2069

Open
devmotion wants to merge 12 commits into
masterfrom
dmw/skewnormal
Open

Add cdf/ccdf/quantile to SkewNormal via Owen's T#2069
devmotion wants to merge 12 commits into
masterfrom
dmw/skewnormal

Conversation

@devmotion

@devmotion devmotion commented Jun 19, 2026

Copy link
Copy Markdown
Member

Stacked on #2070. This PR is based on #2070 (Roots-based quantile root-finding) and
targets dmw/roots; its diff shows only the SkewNormal feature. The near-zero-root
convergence that SkewNormal needs (its support straddles 0) is provided by #2070's solver,
so this PR no longer touches src/quantilealgs.jl. Merge #2070 first.

Summary

SkewNormal was missing cdf/ccdf (the source had #cdf requires Owen's T function.). StatsFuns 2.2 now provides owens_t, which unblocks the closed form:

cdf(x)  = Φ(z) - 2·T(z, α)
ccdf(x) = Φ̄(z) + 2·T(z, α),   z = (x - ξ)/ω

This PR adds cdf/ccdf, refactors pdf/logpdf for consistency, and adds quantile/cquantile. logcdf/logccdf/invlogcdf/invlogccdf use the generic fallbacks.

  • Requires StatsFuns 2.2; bumps the package version.

Quantile root-finding (handled by #2070)

quantile/cquantile are Newton iterations from the mode via quantile_newton/cquantile_newton. The original version of this PR also patched src/quantilealgs.jl with an absolute-tolerance floor, because the hand-rolled relative-only stopping rule never converges as the root approaches zero (which SkewNormal hits). That fix is now subsumed by #2070, which replaces the hand-rolled loops with Roots.find_zero and so converges for near-zero roots out of the box — this PR therefore drops its own quantilealgs.jl change.

invlogcdf/invlogccdf are intentionally left to the generic quantile(d, exp(lp)) fallback rather than the @quantile_newton macro: the log-space invlogcdf_newton overshoots into the underflowed tail from the approximate mode m_0 and returns NaN for extreme log-probabilities (e.g. p ≲ 1e-8), even with #2070's bracketing fallback, whereas routing through quantile_newton stays robust there.

Tests

  • Added a sn-based reference class test/ref/continuous/skewnormal.R (moments from sn:::sn.cumulants, an independent computation rather than a reimplementation of the Julia formulas), 4 entries to continuous_test.lst, the generated entries to continuous_test.ref.json, and pinned sn + deps in renv.lock.
  • The dedicated testset is reduced to construction, support boundaries, and the α=0 reduction to Normal.

Verification

  • R reference framework: 88332/88332 pass for the SkewNormal entries.
  • Reference values cross-checked: R sn and Julia agree to ~1e-14.
  • quantile/cquantile round-trips and the invlogcdf tail fallback verified against Simplify quantile root-finding to Roots defaults #2070's solver (incl. roots straddling 0).

🤖 Generated with Claude Code

Cedriq1astaken and others added 5 commits June 9, 2026 16:06
- 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-commenter

codecov-commenter commented Jun 19, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 86.89%. Comparing base (6b5aa0d) to head (9632e69).

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.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

devmotion and others added 2 commits June 20, 2026 00:50
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>
@devmotion devmotion changed the base branch from master to dmw/roots June 19, 2026 22:53
Base automatically changed from dmw/roots to master June 25, 2026 13:37
@devmotion devmotion marked this pull request as ready for review June 26, 2026 12:36
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants