Statistics on truncated distributions (`cdf`, `median`, `quantile`) can give plausible but wrong results

[yurivish]

This seems wrong:

julia> using Distributions

julia> td = truncated(Normal(), 8, 9);

julia> cdf(td, 8.1)                         # true 0.5583754014201245 (27.87% too high)
0.7139890769706131

julia> ccdf(td, 8.1)                        # true 0.4416245985798754 (19.16% too low)
0.35699453848530655

julia> median(td)                           # true 8.084888899018166 (0.10% too low)
8.07657100413013

julia> quantile(td, 0.9)                    # true 8.278609037011552 (0.83% too low)
8.209536151601391

The error can be subtle, or at other times quite large:

julia> td = truncated(Normal(), 9, 10);

julia> cdf(td, 9.1)                         # true 0.5998422102634853 (off by 100%)
0.0

julia> quantile(td, 0.9)                    # true 9.24935337636521
10.0

Possible cause (summarized by Sonnet)

Truncated caches lcdf = cdf(d0, lower) and tp = ucdf - lcdf, then computes:

cdf(d::Truncated, x)      = clamp((cdf(d.untruncated, x) - d.lcdf) / d.tp, 0, 1)
quantile(d::Truncated, p) = quantile(d.untruncated, d.lcdf + p * d.tp)

When lower sits far enough into d0's right tail, lcdf is within eps()/2 of 1.0 and rounds to exactly 1.0 as a Float64 — even though the true value is 1 - tp for a tp that's still nonzero (down to ~1e-300, way above Float64's underflow floor). Once that rounding happens:

• cdf(d.untruncated, x) - d.lcdf subtracts two numbers that agree to the last bit → 0.0 for every x inside the window.
• d.lcdf + p * d.tp is 1.0 + (tiny) == 1.0 for every p → quantile calls quantile(d0, 1.0), which is Inf by definition, then that Inf gets clamped to upper (or stays Inf if upper == Inf).

So it's pure floating-point cancellation, not a logic error — the formulas are exactly right when lcdf/tp are computed in exact arithmetic.

Why it's fixable without a design change: the struct already stores loglcdf/logtp, and logcdf/logccdf use logsubexp on those instead of subtracting lcdf, so they're correct on the same object. rand (l. 218–233) already has the fix for this exact issue — it switches to computing in log space via invlogcdf(d0, logaddexp(d.loglcdf, d.logtp - randexp(rng))) once tp <= sqrt(eps()), citing #1548/#1553 (a previously-fixed instance of this same cancellation, but only for rand).

The fix: apply the same log-space trick to the other four functions:

cdf(d, x)      = exp(logcdf(d, x))
ccdf(d, x)     = exp(logccdf(d, x))
quantile(d, p) = invlogcdf(d.untruncated, logaddexp(d.loglcdf, d.logtp + log(p)))
cquantile(d, p) = invlogccdf(d.untruncated, logaddexp(d.loglccdf, d.logtp + log(p)))  # symmetric, needs loglccdf cached
median(d)      = quantile(d, 0.5)  # falls out for free once quantile is fixed

or equivalently cache uccdf/lccdf (the ccdf-based tail values) and pick whichever side (lcdf-based vs ccdf-based) avoids cancellation depending on whether lcdf > 1/2.