From 325d36a1c03800d8ff60522d4d8a99502a07c4da Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 14:26:25 +0300 Subject: [PATCH 1/3] Skip PSIS tail smoothing when a quarter of the tail ties at its minimum --- src/stan/services/pathfinder/psis.hpp | 23 +++++-- .../unit/services/pathfinder/psis_test.cpp | 68 +++++++++++++++++++ 2 files changed, 87 insertions(+), 4 deletions(-) diff --git a/src/stan/services/pathfinder/psis.hpp b/src/stan/services/pathfinder/psis.hpp index b4da2550639..1c58ee01f5b 100644 --- a/src/stan/services/pathfinder/psis.hpp +++ b/src/stan/services/pathfinder/psis.hpp @@ -57,7 +57,8 @@ inline Eigen::Array profile_loglikelihood( * @param[in] min_grid_pts The minimum number of grid points used in the fitting * algorithm. * @return A pair of doubles with the first element `sigma` and the second - * element `k`. + * element `k`. Both are NaN if the first quartile of `x` is not larger than + * the smallest value of `x`. * * @details Here the parameter `k is the negative of `k` in Zhang & Stephens * (2009). @@ -80,6 +81,13 @@ inline std::pair gpdfit(const EigArray& x, const double x_1st_qt = x_ref.coeff( static_cast(std::floor(static_cast(N) / 4.0 + 0.5)) - 1l); + // At least a quarter of x equals the minimum. The sample is then far from a + // generalized Pareto distribution, and with a minimum of 0 the grid for + // theta below divides by 0. + if (!(x_1st_qt > x_ref.coeff(0))) { + constexpr double nan = std::numeric_limits::quiet_NaN(); + return {nan, nan}; + } array_vec_t theta = 1.0 / x_ref.coeff(N - 1) + (1.0 - (M / (linspaced_arr - 0.5)).sqrt()) / (prior * x_1st_qt); @@ -121,7 +129,8 @@ inline auto qgpd(const EigArray& p, const double k, const double sigma) { * @param[in] cutoff * @return A pair containing: * `first`: Eigen Array same size as `x` containing the logs of the - * order statistics of the generalized pareto distribution. + * order statistics of the generalized pareto distribution, or `x` itself + * if the shape parameter estimate is not finite. * `second`: scalar shape parameter estimate. */ template @@ -129,7 +138,7 @@ inline auto psis_smooth_tail(const EigArray& x, const double cutoff) { const double exp_cutoff = std::exp(cutoff); const auto fit = gpdfit(x.array().exp() - exp_cutoff); const double k = fit.second; - if (!std::isinf(k)) { + if (std::isfinite(k)) { const Eigen::Index x_size = x.size(); const double sigma = fit.first; auto p @@ -262,7 +271,13 @@ inline Eigen::Array psis_weights( for (Eigen::Index i = 0; i < idx_size; ++i) { llr_weights.coeffRef(idx.coeff(i)) = smoothed.first.coeff(i); } - if (smoothed.second > 0.7) { + if (!std::isfinite(smoothed.second)) { + logger.warn( + "In PSIS Weight Calculation: At least a quarter of the importance " + "ratios in the tail are equal to the smallest one, which is too " + "many for estimating the generalized pareto values. Returning " + "non-pareto smoothed weights."); + } else if (smoothed.second > 0.7) { std::stringstream s; s << "Pareto k value (" << std::setprecision(2) << smoothed.second << ") is greater than 0.7. Importance resampling was not able to " diff --git a/src/test/unit/services/pathfinder/psis_test.cpp b/src/test/unit/services/pathfinder/psis_test.cpp index 9a77d4fb868..82262d9e0e9 100644 --- a/src/test/unit/services/pathfinder/psis_test.cpp +++ b/src/test/unit/services/pathfinder/psis_test.cpp @@ -183,3 +183,71 @@ TEST(ServicesPSIS, max_n_elements) { EXPECT_EQ(sorted_idx(i), sorted_result_pos(i)); } } + +TEST(ServicesPSIS, gpdfit_tied_first_quartile) { + Eigen::Array x(20); + x << 0.00231135747917145, 0.00433831801177895, 0.0108541508266367, + 0.0146361066006147, 0.016979809437058, 0.0175260143161184, + 0.0185788657942587, 0.0221000771645818, 0.0246924617215394, + 0.0274311810928377, 0.0492474466172383, 0.0623614301581945, + 0.0659810585631264, 0.0796802961280512, 0.146305816395337, + 0.184840740340265, 0.20479110080993, 0.257155687804798, 0.414031661251632, + 0.95242504763768; + // 4 of 20 values at the minimum: the first quartile (the 5th value) is + // larger, so the fit is finite + x.head(4).setZero(); + auto fit = stan::services::psis::internal::gpdfit(x); + EXPECT_TRUE(std::isfinite(std::get<0>(fit))); + EXPECT_TRUE(std::isfinite(std::get<1>(fit))); + // 5 of 20 values at the minimum, which is 0 + x.head(5).setZero(); + fit = stan::services::psis::internal::gpdfit(x); + EXPECT_TRUE(std::isnan(std::get<0>(fit))); + EXPECT_TRUE(std::isnan(std::get<1>(fit))); + // 5 of 20 values at a positive minimum + x.head(5).setConstant(0.001); + fit = stan::services::psis::internal::gpdfit(x); + EXPECT_TRUE(std::isnan(std::get<0>(fit))); + EXPECT_TRUE(std::isnan(std::get<1>(fit))); +} + +TEST(ServicesPSIS, psis_weights_tied_cutoff) { + // The tail has 20 values: 15 spread from -3.5 to 0 and 5 of the 15 values + // at -5. The cutoff is also -5, so 5 tail values have a zero excess. + Eigen::Array lr(100); + for (Eigen::Index i = 0; i < 70; ++i) { + lr(i) = -20.0 - 0.1 * i; + } + lr.segment(70, 15).setConstant(-5.0); + for (Eigen::Index i = 85; i < 100; ++i) { + lr(i) = -0.25 * (i - 85); + } + std::unique_ptr log(new std::stringstream{}); + stan::test::test_logger logger(std::move(log)); + auto w = stan::services::psis::psis_weights(lr, 20, logger); + for (Eigen::Index i = 0; i < lr.size(); ++i) { + EXPECT_FLOAT_EQ(w(i), std::exp(lr(i))) << "for weight " << i; + } + EXPECT_NE(logger.log_->str().find("Returning non-pareto smoothed weights"), + std::string::npos); +} + +TEST(ServicesPSIS, psis_weights_underflow_in_tail) { + // Ratios more than about 745 below the largest one underflow to 0 in the + // generalized Pareto fit. Here 5 of the 20 tail values and the cutoff do. + Eigen::Array lr(100); + for (Eigen::Index i = 0; i < 85; ++i) { + lr(i) = -1000.0 - i; + } + for (Eigen::Index i = 85; i < 100; ++i) { + lr(i) = -0.25 * (i - 85); + } + std::unique_ptr log(new std::stringstream{}); + stan::test::test_logger logger(std::move(log)); + auto w = stan::services::psis::psis_weights(lr, 20, logger); + for (Eigen::Index i = 0; i < lr.size(); ++i) { + EXPECT_FLOAT_EQ(w(i), std::exp(lr(i))) << "for weight " << i; + } + EXPECT_NE(logger.log_->str().find("Returning non-pareto smoothed weights"), + std::string::npos); +} From 51712709d200075b18e5ae6bb260bfb89b319a37 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 16:31:30 +0300 Subject: [PATCH 2/3] doc: fix capitalization of Pareto (no code change) --- src/stan/services/pathfinder/psis.hpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/stan/services/pathfinder/psis.hpp b/src/stan/services/pathfinder/psis.hpp index 1c58ee01f5b..00fa05ee0b7 100644 --- a/src/stan/services/pathfinder/psis.hpp +++ b/src/stan/services/pathfinder/psis.hpp @@ -14,16 +14,16 @@ namespace psis { namespace internal { /** - * Compute log joint likelihood parameter estimates from generalized pareto + * Compute log joint likelihood parameter estimates from generalized Pareto * distribution and the samples the parameters were estimated from. * @tparam EigArray1 An Eigen type inheriting from `ArrayBase` with dynamic * compile time rows and 1 compile time column. * @tparam EigArray2 An Eigen type inheriting from `ArrayBase` with dynamic * compile time rows and 1 compile time column. - * @param[in] theta Estimates from generalized pareto distribution estimation + * @param[in] theta Estimates from generalized Pareto distribution estimation * @param[in] x The sample that the parameters were estimated from. * @return Array of the joint log likelihood of parameter estimates from - * generalized pareto distribution and the samples the parameters were estimated + * generalized Pareto distribution and the samples the parameters were estimated * from. */ template @@ -41,10 +41,10 @@ inline Eigen::Array profile_loglikelihood( } /** - * Estimate parameters of the Generalized Pareto distribution + * Estimate parameters of the generalized Pareto distribution * * Given a sample `x`, Estimate the parameters `k` and $sigma$ of - * the Generalized Pareto Distribution (GPD), assuming the location parameter is + * the generalized Pareto distribution (GPD), assuming the location parameter is * 0. By default the fit uses a prior for `k`, which will stabilize * estimates for very small sample sizes (and low effective sample sizes in the * case of MCMC samples). The weakly informative prior is a Gaussian prior @@ -105,7 +105,7 @@ inline std::pair gpdfit(const EigArray& x, } /** - * Inverse CDF of generalized pareto distribution + * Inverse CDF of generalized Pareto distribution * (assuming location parameter is 0) * * @tparam EigArray An Eigen type inheriting from `ArrayBase` with dynamic @@ -129,7 +129,7 @@ inline auto qgpd(const EigArray& p, const double k, const double sigma) { * @param[in] cutoff * @return A pair containing: * `first`: Eigen Array same size as `x` containing the logs of the - * order statistics of the generalized pareto distribution, or `x` itself + * order statistics of the generalized Pareto distribution, or `x` itself * if the shape parameter estimate is not finite. * `second`: scalar shape parameter estimate. */ @@ -262,8 +262,8 @@ inline Eigen::Array psis_weights( std::string("In PSIS Weight Calculation: Difference " "between the tails is ") + std::to_string(eps_diff) + - " which is too small for estimating the generalized pareto values." - " Returning non-pareto smoothed weights."); + " which is too small for estimating the generalized Pareto values." + " Returning non-Pareto smoothed weights."); } else { auto smoothed = internal::psis_smooth_tail(lw_tail, cutoff); auto idx = max_n.second.tail(tail_len); @@ -275,8 +275,8 @@ inline Eigen::Array psis_weights( logger.warn( "In PSIS Weight Calculation: At least a quarter of the importance " "ratios in the tail are equal to the smallest one, which is too " - "many for estimating the generalized pareto values. Returning " - "non-pareto smoothed weights."); + "many for estimating the generalized Pareto values. Returning " + "non-Pareto smoothed weights."); } else if (smoothed.second > 0.7) { std::stringstream s; s << "Pareto k value (" << std::setprecision(2) << smoothed.second From daa080dc479d8ea9c701c854c326231a5c0082e4 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 19 Oct 2026 17:19:59 +0300 Subject: [PATCH 3/3] Match the capitalized PSIS warning in the tests --- src/test/unit/services/pathfinder/psis_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/test/unit/services/pathfinder/psis_test.cpp b/src/test/unit/services/pathfinder/psis_test.cpp index 82262d9e0e9..e6ccd544e16 100644 --- a/src/test/unit/services/pathfinder/psis_test.cpp +++ b/src/test/unit/services/pathfinder/psis_test.cpp @@ -228,7 +228,7 @@ TEST(ServicesPSIS, psis_weights_tied_cutoff) { for (Eigen::Index i = 0; i < lr.size(); ++i) { EXPECT_FLOAT_EQ(w(i), std::exp(lr(i))) << "for weight " << i; } - EXPECT_NE(logger.log_->str().find("Returning non-pareto smoothed weights"), + EXPECT_NE(logger.log_->str().find("Returning non-Pareto smoothed weights"), std::string::npos); } @@ -248,6 +248,6 @@ TEST(ServicesPSIS, psis_weights_underflow_in_tail) { for (Eigen::Index i = 0; i < lr.size(); ++i) { EXPECT_FLOAT_EQ(w(i), std::exp(lr(i))) << "for weight " << i; } - EXPECT_NE(logger.log_->str().find("Returning non-pareto smoothed weights"), + EXPECT_NE(logger.log_->str().find("Returning non-Pareto smoothed weights"), std::string::npos); }