[avehtari]

This PR was assisted by Claude. The reference code in R was written by human.

Summary

Zhang and Stephens (2009) method requires that first quartile is bigger than the smallest value. psis_smooth_tail() passes the excesses exp(lw_tail) - exp(cutoff) to gpdfit(). When a quarter of them are 0, the first quartile x_1st_qt is 0, and the grid theta divides by it and becomes -Inf. The product -theta * x is then NaN for the zero excesses, so k and sigma are NaN. The test !std::isinf(k) lets NaN through, the smoothed tail is NaN, and (llr_weights.array() < 0.0).select(llr_weights, 0.0) maps NaN to 0. Each tail draw gets the weight exp(0) = 1. The k > 0.7 warning does not fire for NaN.

Zero excesses occur in these cases:

• tail ratios tie with the cutoff, for example duplicate draws;
• the cutoff and some tail ratios are more than about 745 below the largest ratio, so exp() underflows to 0 for them;
• many draws have a log density of -Inf, which single.hpp sets when the log density throws a std::domain_error.

Example (100 draws, tail length 20): 70 log ratios from -20 to -26.9, 15 at -5, and 15 from -3.5 to 0. The cutoff and 5 tail values are -5. develop returns 1 for all 20 tail draws; the correct unsmoothed weights are exp(lr), for example exp(-5) = 0.0067 for the 5 tied draws.

When the tied minimum is positive (a quarter of the tail ties above the cutoff), theta is finite, and gpdfit() fits a sample that is far from a generalized Pareto distribution. posterior::gpdfit() returns NA for both cases, with the same check on the first quartile.

Intended Effect

• gpdfit() returns NaN for sigma and k when the first quartile of x is not larger than its smallest value.
• psis_smooth_tail() smooths only for a finite k. Otherwise it returns the tail unchanged.
• psis_weights() logs a warning when k is not finite and returns the weights without smoothing, as it already does for a constant tail.

This behavior matches the behavior of corresponding function in posterior R package.

How to Verify

./runTests.py src/test/unit/services/pathfinder. The three new tests in psis_test.cpp are:

• gpdfit_tied_first_quartile: 4 of 20 values at the minimum give a finite fit; 5 of 20 at a minimum of 0, and 5 of 20 at a minimum of 0.001, give NaN.
• psis_weights_tied_cutoff: the example above. The weights must equal exp(lr), and the warning must be logged.
• psis_weights_underflow_in_tail: 85 log ratios from -1000 to -1084 and 15 from -3.5 to 0, so 5 tail values and the cutoff underflow. The same checks.

With the develop psis.hpp and the new tests, the 3 new tests fail (42 expectations) and the 5 existing tests pass. With this PR all 8 pass. eight_schools_test and normal_glm_test also pass.

Tested with develop 849cec8da and Stan Math e4745c87ff, g++ 14.3.0, Linux x86_64.

Side Effects

The results do not change when the first quartile of the excesses is larger than their minimum. The existing PSIS tests and the Pathfinder service tests pass without changes. A tail with a quarter of its values tied at a positive minimum is no longer smoothed.

Documentation

• The doc comments of gpdfit() and psis_smooth_tail() describe the new return values.
• The second commit fixes the capitalization of Pareto in docs, warnings and comments (no code changes)

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

Aalto University

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