Skip PSIS tail smoothing when a quarter of the tail ties at its minimum

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 <typename EigArray1, typename EigArray2>
@@ -41,10 +41,10 @@ inline Eigen::Array<double, Eigen::Dynamic, 1> 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
@@ -57,7 +57,8 @@ inline Eigen::Array<double, Eigen::Dynamic, 1> 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<double, double> gpdfit(const EigArray& x,
   const double x_1st_qt = x_ref.coeff(
       static_cast<Eigen::Index>(std::floor(static_cast<double>(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<double>::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);
@@ -97,7 +105,7 @@ inline std::pair<double, double> 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
@@ -121,15 +129,16 @@ 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 <typename EigArray>
 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
@@ -253,16 +262,22 @@ inline Eigen::Array<double, Eigen::Dynamic, 1> 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);
       const Eigen::Index idx_size = idx.size();
       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 "

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<double, -1, 1> 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<double, -1, 1> 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<std::stringstream> log(new std::stringstream{});
+  stan::test::test_logger<std::stringstream> 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<double, -1, 1> 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<std::stringstream> log(new std::stringstream{});
+  stan::test::test_logger<std::stringstream> 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);
+}