Centralize numerical quantile inversion for internal distributions

[lrnv]

Problem

Several internal univariate distributions reimplement numerical quantiles with nearly identical ad-hoc root finding. The duplication currently covers multiple subtly different versions:

• finite-support bisection in ClaytonWilliamsonDistribution;
• positive unbounded-support bracket expansion in WilliamsonBetaProduct and WilliamsonFromFrailty;
• a separate 64-step loop in PowerTiltedFrailty;
• unit-interval inversion, sometimes through logcdf, for Distortion;
• specialized segmented/cache-aware inversion in LiouvilleConditionalRadial.

This duplicates endpoint validation, support handling, bracket growth, stopping rules, floating-point type selection, and treatment of inaccurate CDF values near boundaries. The recent Liouville failures show that these details are easy to get inconsistently wrong.

Related organizational audit: #420.

Proposed interface

Provide one internal, support-aware numerical quantile algorithm. Each distribution keeps an explicit one-line Distributions.quantile method, avoiding any broad extension on Distributions.UnivariateDistribution:

Distributions.quantile(d::MyDistribution, p::Real) =
    _quantile_from_cdf(quantile_strategy(typeof(d)), d, p)

The distribution instance already carries its model parameters. The shared interface must receive only numerical-strategy information.

Internal strategy trait

Use a small internal trait hierarchy to select the numerical objective consistently for each distribution type:

abstract type QuantileStrategy end
struct CDFQuantile <: QuantileStrategy end
struct LogCDFQuantile <: QuantileStrategy end

quantile_strategy(::Type) = CDFQuantile()
quantile_strategy(::Type{<:SomeDistortion}) = LogCDFQuantile()

The strategy then dispatches into the common engine:

_quantile_from_cdf(::CDFQuantile, d, p) = ...
_quantile_from_cdf(::LogCDFQuantile, d, p) = ...

Additional orthogonal traits, for example bounded versus expanding brackets, should be introduced only if the audit demonstrates that support information from minimum(d) and maximum(d) is insufficient. Prefer deriving the bracket behavior directly from the support.

No macro is needed. Do not redefine quantile for all Distributions.UnivariateDistribution objects: both the function and that abstract type belong to Distributions.jl, and such a broad method could interfere with its own fallbacks or downstream packages.

Keep the trait and helper internal under YAGNI.

Required semantics

The common implementation should compute the generalized quantile

inf{x : F(x) >= p}

rather than assume a smooth strictly increasing CDF. It should centralize:

• validation of p in [0, 1];
• exact p = 0 and p = 1 support endpoints;
• finite, semi-infinite, and possibly fully unbounded supports;
• robust bracket expansion;
• cdf versus logcdf objectives where underflow matters;
• type-aware stopping rules and return types;
• atoms, plateaus, and finite endpoint masses;
• NaN/non-finite objective diagnostics.

Closed-form quantiles and genuinely structured algorithms, such as the cached segmented inversion of LiouvilleConditionalRadial, should remain specialized.

Audit scope

Inventory every Distributions.quantile method and every root search used as a quantile. Classify each as:

1. closed form;
2. common numerical inversion candidate;
3. specialized numerical inversion with a documented reason.

Initial common candidates include WilliamsonFromFrailty, PowerTiltedFrailty, ClaytonWilliamsonDistribution, WilliamsonBetaProduct, generic Williamson inverse distributions, Distortion fallbacks, and BernsteinDistortion through _unit_quantile.

Acceptance criteria

• One tested numerical inversion engine covers the ordinary finite/unit/positive-unbounded cases.
• A new distribution opts in through an explicit one-line quantile method.
• An internal trait selects cdf versus logcdf without repeating keyword choices at call sites.
• Support and bracketing are inferred from the distribution whenever possible.
• Closed-form and cache-aware quantiles keep their specialized dispatch.
• Tests cover bounded and unbounded support, Float32/Float64 where practical, p = 0/1, atoms, plateaus, endpoint mass, extreme probabilities, and failure to establish a valid bracket.
• Existing duplicated implementations are migrated only when behavior is demonstrably equivalent.