Remove the offset at a tied Pareto tail cutoff

[avehtari]

This PR was assisted by Claude

Summary

Offset at a tied Pareto tail cutoff (added in d9f7627) is my mistake.

When the largest non-tail value ties with the smallest tail values, ps_tail() and pareto_pit() subtract .Machine$double.eps from the cutoff. This offset causes three problems:

1. pareto_smooth() can return values below the cutoff. For non-negative weights with many exact zeros, the result contains negative weights of about -2e-16, and resample_draws() then fails with Assertion on 'weights' failed: Element <i> is not >= 0.
2. pareto_pit() puts an observation at a tied cutoff, or less than eps below it, in the GPD tail. For example, with 350 of 400 draws at 0, pareto_pit() gives 0.85 for y = -1e-300, which is below all draws.
3. The offset is absolute, so the fit depends on the scale of the data. It has no effect when |cutoff| >= 2, and it changes k-hat when the tail excesses are small. For natural-scale weights that span many orders of magnitude, k-hat differs from the fit to the same weights on the log scale.

This PR removes the offset. Tail values equal to the cutoff then give zero excesses, which the GPD fit and the GPD quantile and distribution functions handle correctly.

Why zero excesses are valid

• The GPD support starts at mu and includes it: the density at mu is 1 / sigma, finite and positive. A sample value at mu has a finite log-likelihood.
• gpdfit() uses the excesses in log1p(-theta * x), which is 0 for x = 0, and in 1 / x[N] and 1 / xstar. x[N] > 0 because the tail is not constant.
• The degenerate case is xstar == x[1] (at least a quarter of the excesses equal the minimum). gpdfit() already returns k = NA for it (be503f2). The offset did not prevent this case: all tied values got the same excess eps.
• qgeneralized_pareto() returns mu plus a value >= 0, so smoothed values are never below the cutoff. No clamp is necessary.
• The code already ran with zero excesses whenever |cutoff| >= 2, because cutoff - eps == cutoff there.

How the negative weights occur in practice

cmdstanr process_init_approx() smooths exp(lp__ - lp_approx__ - max(...)) for Pathfinder inits. When most Pathfinder paths are poor, most of these values underflow to exactly 0, and the cutoff and some tail values are both 0. The negative weights then go to resample_draws(), and $sample(init = pathfinder_fit) fails. This was observed as a deterministic failure in a cluster run with 10 paths of 40 draws.

Minimal example (no Stan):

library(posterior)
w <- c(rep(0, 350), exp(-seq(0, 60, length.out = 50)))
ws <- pareto_smooth(w, tail = "right", r_eff = 1, return_k = TRUE, verbose = FALSE)
sum(ws$x < 0)          # 7 on master, 0 with this PR
ws$diagnostics$khat    # 8.86 on master, 18.81 with this PR
ps_tail(log(w), ndraws_tail = 60, are_log_weights = TRUE)$k  # log scale: 18.81 on both
x <- as_draws_df(data.frame(theta = seq_along(w)))
resample_draws(x, ndraws = 4, weights = ws$x, method = "simple_no_replace")
#> master: Error: Assertion on 'weights' failed: Element 341 is not >= 0.

pareto_pit() example:

draws <- c(rep(0, 350), exp(-seq(0, 5, length.out = 50)))
x0 <- as_draws_matrix(matrix(draws, ncol = 1))
x5 <- as_draws_matrix(matrix(draws + 5, ncol = 1))
pareto_pit(x0, -1e-300)          # master 0.85, this PR 2.5e-07 (below all draws)
set.seed(1); pareto_pit(x0, 0)   # master 0.85 (GPD tail), this PR 0.2323 (randomized raw PIT)
set.seed(1); pareto_pit(x5, 5)   # 0.2323 on master and with this PR, because 5 - eps == 5

Scale dependence on master, x = s * c(rep(0, 350), exp(rnorm(50, 0, 2))) with set.seed(3) and 60 tail draws: k-hat is 1.6771 for s = 1 and s = 1e-8, 1.6694 for s = 1e-12, 1.3330 for s = 1e-14, and 0.1066 for s = 1e-16. With this PR, k-hat is 1.6771 for all these s.

Change

• ps_tail() (R/pareto_smooth.R): remove the offset. Update the comment on the cutoff.
• pareto_pit.draws_matrix() (R/pit.R): remove the offset in the right and the left tail.

Tests

• tests/testthat/test-pareto_smooth.R:
  • ps_tail does not smooth tail values below the cutoff: right tail and left tail, ties at 0. Checks the sign and that the order of the draws does not change.
  • pareto_smooth does not return negative weights for ties at 0: checks the sign and that resample_draws() accepts the weights.
• tests/testthat/test-pareto_pit.R:
  • pareto_pit uses the raw PIT at a tail cutoff that ties with tail draws: right and left tail. Checks that draws tied at 0 and draws tied at 5 give the same PIT for an observation at the tie, and that an observation outside all draws gets the smallest or largest PIT.

On unpatched master, the new tests give 9 failures: 5 in test-pareto_smooth.R (one of them the resample_draws() assertion error) and 4 in test-pareto_pit.R. With this PR, all tests pass.

Effect on other input

ps_tail() on master and with this PR on 1200 random inputs (100, 400 and 1000 draws):

input	identical	max abs. difference
normal draws, 400 runs	400	0
log-normal draws, 400 runs	400	0
log weights, 200 runs	200	0
log weights with -Inf, 200 runs	200	0
rounded normal draws (ties at non-zero cutoffs), 400 runs	60	5.8e-15
natural-scale weights with ties at 0, 400 runs	244	0.98

• In the last row, master gave 1180 negative values, and this PR gave none. The large differences come from the change of k-hat. With this PR, the natural-scale result equals the log-scale result for the same weights: in 200 cases, the smoothed values agree to 1.9e-14, k-hat agrees to 1.6e-13, and k-hat is NA in the same 57 cases.
• In a Pathfinder-like design (350 log weights below -1000, 50 from N(0, 20), 200 seeds), master gave k-hat = NA in 90 seeds; this PR gave NA in none. The offset made tail values below about 2e-16 equal to the tied zeros, which triggered the xstar == x[1] case.

Tested versions

• posterior: this branch, based on upstream master 6dd1e12383.
• R 4.6.1, testthat 3.3.2, Linux x86_64.
• devtools::test(): 2232 tests, 0 failures, 0 errors, 4 skips (gbm is not installed).

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/)

[github-actions]

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

• ✔️as_draws_array: 155ms -> 154ms [-1.13%, +0.46%]
• ✔️as_draws_df: 77.3ms -> 77.5ms [-0.9%, +1.45%]
• ✔️as_draws_list: 172ms -> 172ms [-0.46%, +1.35%]
• ✔️as_draws_matrix: 25.7ms -> 25.8ms [-0.95%, +1.88%]
• ✔️as_draws_rvars: 126ms -> 126ms [-0.9%, +1.54%]
• ✔️summarise_draws_100_variables: 762ms -> 760ms [-0.74%, +0.27%]
• ✔️summarise_draws_10_variables: 85.3ms -> 85.5ms [-0.58%, +1.01%]
  Further explanation regarding interpretation and methodology can be found in the documentation.

[florence-bockting]

Thank you. Looks good. I just added a bullet point to NEWS.md.

[avehtari]

@florence-bockting your NEWS commit caused a conflict, can you resolve it?

[florence-bockting]

@florence-bockting your NEWS commit caused a conflict, can you resolve it?

Done