Return the support limit from ssd_q*() at the extreme quantiles

[joethorley]

Closes #195

Summary

.qd() short-circuited p = 0 and p = 1 to 0 and Inf before dispatching to the distribution's own quantile function. That assumes every distribution has support unbounded above and approaching zero below, which is true for most of the set but not for invpareto, bounded above at scale:

ssd_qinvpareto(1, shape = 1, scale = 5)             #> Inf, should be 5
ssdtools:::qinvpareto_ssd(1, shape = 1, scale = 5)  #> 5

So the exported function discarded the correct answer its own backend had already computed, and did not round trip through ssd_pinvpareto() at the upper support limit.

Investigating why the shortcut was there turned up the real cause. root() in R/internal.R uses uniroot(lower = 0, upper = 1, extendInt = "upX"), which extends the bracket rather than diverging, so the three quantile functions that solve numerically saturate at the endpoints instead of returning infinity. The shortcut masked that, at the cost of overriding the distributions with genuinely finite support.

Measured on the live .qd() dispatch targets before this change:

dispatch target	p = 0	p = 1	correct
burrIII3, gamma, gompertz, lnorm, weibull	0	Inf	yes
gumbel, logis (both .lgt)	-Inf	Inf	yes
invpareto	0	scale	yes, finite upper
logis_logis (.lgt)	-1073.742	41.95	no, should be ∓Inf
lnorm_lnorm	0	10486.75	no, should be Inf
multi	errors	5243.9	no, should be 0 and Inf

Eight of eleven targets already returned their own correct endpoints unaided.

What changed

endpoints() in R/internal.R root finds the interior probabilities and substitutes the exact quantiles at p = 0 and p = 1. It is used by qlogis_logis_ssd() (-Inf, Inf on the log scale), qlnorm_lnorm_ssd() (0, Inf) and qmulti_list().

qmulti_list() derives its limits from the components rather than hard coding (0, Inf): the mixture's support is the union of the positive-weight components' supports, and each component's q*_ssd() now reports its own limits exactly, so lower and upper are the min and max of those at p = 0 and p = 1. A mixture weighted entirely on a bounded component returns its finite limit, and any future bounded distribution inherits correct model-averaged endpoints without touching this code.

The p == 0 / p == 1 block is then removed from .qd(), so qdist() always defers to the distribution. .lgt is no longer read in .qd() or .qdist() and was dropped from both; qdist() still uses it to back transform.

pinvpareto_ssd() is clamped to [0, 1], the p-side counterpart of the same bounded-support contract. It returned (q / scale)^shape unclamped, so the CDF exceeded 1 above scale and was negative for negative q with odd shape: ssd_pinvpareto(10, shape = 3, scale = 5) was 8, and ssd_hp() on the ccme_boron inverse Pareto fit (fitted scale 75.3) reported 117.5% of species affected at a concentration of 100.

test_dist() gains lower and upper arguments, because it asserted q(0) == 0 and q(1) == Inf for every distribution, which is the same assumption in test form. Its CDF monotonicity probe now sits strictly inside the support rather than at q = 1, which for invpareto with default scale = 1 is the flat region above the limit. Its ssd_qmulti() endpoint assertions now supply lnorm.weight = 1; they had relied on the shortcut returning before normalize_weights() could raise "at least one distribution must have a positive weight". A new test-multi.R block checks that a single-component mixture reproduces that component's endpoints for every distribution in ssd_dists_all(), and that a mixture on a bounded component returns its finite limit.

This supersedes the direction suggested in #195, which proposed a support flag in dist_data. That would change a user-visible dataset and would leave the three numerical quantile functions still wrong when called directly, so it treats the symptom rather than the cause.

Behaviour changes

Beyond the fixes above, three exported behaviours change at the endpoints, all arguably corrections:

• ssd_hc() accepts proportion of 0 and 1 (chk_range() is inclusive), so for a model average these now flow through qmulti_list()'s derived limits. For the current distribution set they are still 0 and Inf.
• With the shortcut gone, invalid parameters at p = 0 or 1 return NaN rather than 0/Inf, consistent with every other p (e.g. ssd_qlnorm(0, sdlog = -1)).
• ssd_qmulti(0) and ssd_qmulti(1) with no positive weight now error like every other p, instead of returning 0/Inf before validation.

Verification

R CMD check reports 0 errors, 0 warnings and 0 notes. The full test suite is 1370 passing with 0 failures and 0 warnings. air format --check . is clean.

Tests were written first and watched fail: 6 failures on the pre-fix code for the endpoint change (reporting Inf where scale was expected and the saturated -1073.7, 42.0 and 10486.8 where infinite endpoints were expected), and 4 for the clamp and derived-support change (8, 1.2 2.0, -0.008 from the unclamped CDF, Inf from the hard-coded mixture limit).

Endpoints across every exported quantile function after the change:

ssd_qburrIII3        q(0) = 0    q(1) = Inf
ssd_qgamma           q(0) = 0    q(1) = Inf
ssd_qgompertz        q(0) = 0    q(1) = Inf
ssd_qinvpareto       q(0) = 0    q(1) = 1     # = scale, was Inf
ssd_qlgumbel         q(0) = 0    q(1) = Inf
ssd_qllogis          q(0) = 0    q(1) = Inf
ssd_qllogis_llogis   q(0) = 0    q(1) = Inf
ssd_qlnorm           q(0) = 0    q(1) = Inf
ssd_qlnorm_lnorm     q(0) = 0    q(1) = Inf
ssd_qweibull         q(0) = 0    q(1) = Inf
ssd_qmulti           q(0) = 0    q(1) = Inf

And the round trip that was broken:

ssd_pinvpareto(5, shape = 1, scale = 5)  #> 1
ssd_qinvpareto(1, shape = 1, scale = 5)  #> 5

Notes

Branched from main and targeted at dev per the fork convention. ltriangle (#169) has bounded support on both sides and inherits this fix when that branch merges; its test_dist("ltriangle") call will then need lower = exp(-3), upper = exp(3), so the merge order is this PR to dev, dev into add-triangle-distribution, then that one-line change there.

The pre-existing cost of rebuilding the mixture skeleton on every scalar .qd() call, which the review of this PR measured, is filed separately as #197.