Improve floating point accuracy

[avehtari]

Summary

Herbie https://herbie.uwplse.org/ checks floating point computations, detects inaccurate expressions and finds more accurate replacements. I used Herbie with Sol to find these fixes. I have checked that all of them make sense.

• gdp.R: Stable upper-tail and log-probability handling

• gdp.R: Treat only an exactly zero GPD shape as exponential. Small nonzero shapes now retain their shape correction and finite-support behavior through the stable log1p() formulation.

• pareto_smooth.R: Stable ps_convergence_rate() around k = 0.5 and k = 1

• pareto_smooth.R: Use expm1() for log-weight differences

• convergence.R: Centered fourth-moment variance calculation

• discrete-summaries.R: Use log1p() in dissent calculation

• New tests

• All new and old tests are passing

There are two new helper functions which are used once or twice, and they could be also inlined, but I think using them improves the readability.

Copyright and Licensing

By submitting this pull request, the copyright holder is agreeing to
license the submitted work under the following licenses:

• Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

Aki Vehtari

[github-actions]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 700d18d is merged into master:

• 🚀as_draws_array: 168ms -> 166ms [-1.98%, -0.98%]
• ❗🐌as_draws_df: 91.1ms -> 92ms [+0.42%, +1.58%]
• ✔️as_draws_list: 177ms -> 176ms [-0.59%, +0.36%]
• ✔️as_draws_matrix: 28.3ms -> 27.9ms [-3.45%, +0.9%]
• ✔️as_draws_rvars: 136ms -> 136ms [-1.08%, +0.59%]
• ✔️summarise_draws_100_variables: 750ms -> 749ms [-0.48%, +0.24%]
• ✔️summarise_draws_10_variables: 83ms -> 83ms [-0.36%, +0.38%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[n-kall]

Nice! I added one comment and one suggestion for inline comment.
Do you have the results from herbie on how much these suggestions improve things?

[avehtari]

Here is the Sol generated summary of the effect of the fixes

Expected benefit by fix

1. GPD tail probabilities and quantiles — high benefit

Affected code:

• qgeneralized_pareto()
• pgeneralized_pareto()

in posterior/R/gpd.R.

Herbie improved the upper-tail log-probability quantile kernel from 53% to 100%.

The original implementation converted a log probability to ordinary scale and then subtracted it
from one. For an upper-tail log probability of -100:

original quantile: Inf
correct quantile:  2425825972.048951

The CDF had the corresponding problem:

exponential GPD at q = 100

original upper probability:     0
correct upper probability:      3.720075976020836e-44
original log upper probability: -Inf
correct log upper probability:  -100

Herbie scored the upper-CDF kernel at 100% even before rewriting. This demonstrates a sampling
limitation: random points rarely hit the narrow region where a representable survival probability is
lost after the lower CDF rounds to one.

The new implementation carries a log-survival probability through the calculation and only converts
to the requested output form at the end.

Expected benefits:

• Preserves extreme but representable upper-tail probabilities.
• Returns finite GPD quantiles for small upper-tail probabilities.
• Makes lower.tail and log.p behavior consistent with standard R distribution functions.
• Improves downstream Pareto-PIT calculations in extreme tails.

2. Pareto convergence rate — high benefit near transition points

Affected code: ps_convergence_rate() in posterior/R/pareto_smooth.R.

Herbie results:

general domain:     99% → 99%
near k = 0.5:       65% → 99%
near k = 1:         70% → 98%

The general score concealed severe cancellation in narrow neighborhoods.

For the next representable double above k = 0.5, with ndraws = 10:

original internal formula: -0.5925925925925926
original public result:     0
high-precision result:      0.6768166292078592

The exact k = 0.5 branch was also discontinuous from the implemented general expression. The new
version uses algebraically equivalent expm1() forms selected according to the nearest singular
limit.

Expected benefits:

• Produces continuous, nonnegative rates around k = 0.5.
• Avoids catastrophic cancellation as k approaches one.
• Correctly returns one at k = 0 and zero at k = 1.
• Improves diagnostic reliability precisely near the thresholds where convergence-rate values are
  most sensitive.

Away from the transition regions, the new and old formulas agreed within approximately
7e-15 in validation sweeps.

3. Pareto smoothing tail differences — moderate to high benefit

Affected code: ps_tail() in posterior/R/pareto_smooth.R.

Herbie improved the focused log-weight exceedance kernel from 7% to 100%.

The original implementation exponentiated shifted log weights and then subtracted the exponentiated
cutoff. For nearby values around -30:

original relative error: 2.29e-4
stable relative error:   4.60e-17

Expected benefits:

• Preserves small positive GPD exceedances near the tail cutoff.
• Improves fitted Pareto shape and scale parameters for clustered log weights.
• Prevents avoidable zero exceedances caused by rounded exponential subtraction.
• Has little effect when tail values are well separated.

4. mcse_sd() fourth-moment calculation — high benefit for large offsets

Affected code: mcse_sd.default() in posterior/R/convergence.R.

Herbie improved the scalar fourth-moment kernel from 60% to 98%. The implemented fix instead
addresses the vector-level source of correlated rounding error.

For 400 symmetric observations near ±1e8:

original fourth-moment difference: 5.40432e16
centered result:                   4e16

original mcse_sd: 0.0575323908112901
centered mcse_sd: 0.0494962056309299
relative error:   16.2%

The new calculation treats squared centered draws as a new variable and computes their variance from
centered deviations.

Expected benefits:

• Avoids cancellation between the fourth moment and squared variance.
• Prevents inaccurate or potentially negative intermediate variance estimates.
• Improves MCSE estimates for large-magnitude draws whose squared magnitudes vary only slightly.
• Leaves ordinary-scale results effectively unchanged.

5. Ordinal dissent logarithm — low but exact benefit

Affected code: dissent.default() in posterior/R/discrete-summaries.R.

Herbie improved the scalar kernel from 8% to 100% by replacing

log2(1 - r)

with the equivalent

log1p(-r) / log(2)

For r = 1e-20:

original result: 0
stable result:  -1.4426950408889633e-20

Expected benefits:

• Retains tiny nonzero category contributions.
• Improves behavior for categories located extremely close to the sample mean.
• Has negligible practical impact for most ordinal samples because larger category contributions
  dominate the total.

[avehtari]

As Herbie didn't consider what happens with probabilities / densities / weights equal to 0 and log of them equal to -Inf, I added edge tests and fixed behavior for edge cases. Full NEWS

• log_sum_exp() shifted by max(max(x), 0) rather than max(x), because
  warnings = FALSE was passed to max(), which has no such argument and so
  treated FALSE as an extra value to maximize over. Sums of log values below
  about -745 underflowed to -Inf: weights() returned Inf for every draw when
  given unnormalized log weights around -1000, and the normalized log weights used
  by pit() and pareto_pit() were affected the same way.
• weight_draws() now errors when the supplied weights carry no probability
  mass, instead of creating an object whose weights() are NaN and whose
  resample_draws() fails with "missing value where TRUE/FALSE needed". The
  message points at log = TRUE, since underflow is the usual cause. This is a
  breaking change.
• weights() now errors rather than returning NaN when every draw has zero
  weight, which can happen after subsetting an object that was valid when built.
  normalize = FALSE still returns zeros.
• pit() and pareto_pit() no longer abort when one variable's weight column
  carries no probability mass. They warn, naming the variable, and return NA for
  it while the remaining variables are computed as usual, matching pareto_khat().
• pgeneralized_pareto() returned -Inf for upper-tail log probabilities that
  are perfectly representable, for example log P(X > 50) with k = 0, because it
  computed the lower tail and then 1 - p. The survival function is now carried on
  the log scale throughout.
• pgeneralized_pareto() treated any |k| < 1e-15 as k = 0. Small nonzero
  shapes are now handled exactly, down to denormal k.
• qgeneralized_pareto() did not validate p: q(-0.5, k = 0.2) returned
  -0.389 and a positive log probability returned a below-support value, both
  silently. Out-of-range probabilities now return NaN with a warning, as in the
  base R quantile functions.
• pareto_convergence_rate() and ps_convergence_rate() were wrong at two exact
  shape values and imprecise near one of them. A Pareto k of exactly 0 returned 0
  rather than 1, so the rate jumped from 1 just below 0 to 0 at 0 and back to ~1
  just above; k of exactly 0.5 dropped a factor of S/(S-1), giving an error of
  16% at S = 10 and 2.8e-4 at S = 4000; and the surrounding region lost about seven
  digits to cancellation. The rate is now computed in a numerically stable form and
  agrees with the Appendix B expression to 3e-15 across k in (0, 1).
• pareto_convergence_rate() and ps_convergence_rate() read a missing Pareto k
  as 0, reporting a rate of 0, i.e. "non-finite mean", for an unknown shape, and
  errored on a vector containing NA. Missing shapes now propagate as NA.
• mcse_sd() computed the variance of the variance as mean(x^4) - mean(x^2)^2,
  which cancels catastrophically for location-shifted draws; it was 16% off on a
  sample centred at 1e8. It is now computed in two passes.
• dissent() evaluated log2(1 - x), which loses all precision for small x
  and returns exactly 0 below about 1e-16.
• ps_tail() aborted inside gpdfit() with "missing value where TRUE/FALSE
  needed" when the tail contained -Inf log weights equal to the cutoff, because
  the tail-excess helper returned NaN for equal infinities.
• pareto_smooth() and ps_tail() compute tail excesses on the log scale, so
  the generalized Pareto fit no longer loses precision when log weights in the tail
  are close to the cutoff.

[florence-bockting]

Alright, I made the changes to the file structure otherwise looks good to me.

[paul-buerkner]

Thank you all! Feel free to merge when you think it's ready.