Merge pull request #3241 from stan-dev/fix-gamma-lccdf

fix-gamma-lccdf

Author: spinkney

stan/math/prim/prob/gamma_lccdf.hpp

@@ -58,53 +58,55 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
   // The gradients are technically ill-defined, but treated as zero
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (y_vec.val(i) == 0) {
+      // LCCDF(0) = log(P(Y > 0)) = log(1) = 0
       return ops_partials.build(0.0);
     }
   }
 
-  VectorBuilder<is_autodiff_v<T_shape>, T_partials_return, T_shape> gamma_vec(
-      math::size(alpha));
-  VectorBuilder<is_autodiff_v<T_shape>, T_partials_return, T_shape> digamma_vec(
-      math::size(alpha));
-
-  if constexpr (is_autodiff_v<T_shape>) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
-      const T_partials_return alpha_dbl = alpha_vec.val(i);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
-    }
-  }
-
   for (size_t n = 0; n < N; n++) {
     // Explicit results for extreme values
     // The gradients are technically ill-defined, but treated as zero
     if (y_vec.val(n) == INFTY) {
+      // LCCDF(∞) = log(P(Y > ∞)) = log(0) = -∞
       return ops_partials.build(negative_infinity());
     }
 
     const T_partials_return y_dbl = y_vec.val(n);
     const T_partials_return alpha_dbl = alpha_vec.val(n);
     const T_partials_return beta_dbl = beta_vec.val(n);
-
-    const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_dbl);
-
-    P += log(Pn);
-
-    if constexpr (is_autodiff_v<T_y>) {
-      partials<0>(ops_partials)[n] -= beta_dbl * exp(-beta_dbl * y_dbl)
-                                      * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                      / tgamma(alpha_dbl) / Pn;
+    const T_partials_return beta_y_dbl = beta_dbl * y_dbl;
+
+    // Qn = 1 - Pn
+    const T_partials_return Qn = gamma_q(alpha_dbl, beta_y_dbl);
+    const T_partials_return log_Qn = log(Qn);
+
+    P += log_Qn;
+
+    if constexpr (is_any_autodiff_v<T_y, T_inv_scale>) {
+      const T_partials_return log_y_dbl = log(y_dbl);
+      const T_partials_return log_beta_dbl = log(beta_dbl);
+      const T_partials_return log_pdf
+          = alpha_dbl * log_beta_dbl - lgamma(alpha_dbl)
+            + (alpha_dbl - 1.0) * log_y_dbl - beta_y_dbl;
+      const T_partials_return common_term = exp(log_pdf - log_Qn);
+
+      if constexpr (is_autodiff_v<T_y>) {
+        // d/dy log(1-F(y)) = -f(y)/(1-F(y))
+        partials<0>(ops_partials)[n] -= common_term;
+      }
+      if constexpr (is_autodiff_v<T_inv_scale>) {
+        // d/dbeta log(1-F(y)) = -y*f(y)/(beta*(1-F(y)))
+        partials<2>(ops_partials)[n] -= y_dbl / beta_dbl * common_term;
+      }
     }
+
     if constexpr (is_autodiff_v<T_shape>) {
+      const T_partials_return digamma_val = digamma(alpha_dbl);
+      const T_partials_return gamma_val = tgamma(alpha_dbl);
+      // d/dalpha log(1-F(y)) = grad_upper_inc_gamma / (1-F(y))
       partials<1>(ops_partials)[n]
-          += grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
-                                digamma_vec[n])
-             / Pn;
-    }
-    if constexpr (is_autodiff_v<T_inv_scale>) {
-      partials<2>(ops_partials)[n] -= y_dbl * exp(-beta_dbl * y_dbl)
-                                      * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                      / tgamma(alpha_dbl) / Pn;
+          += grad_reg_inc_gamma(alpha_dbl, beta_y_dbl, gamma_val, digamma_val)
+             / Qn;
     }
   }
   return ops_partials.build(P);

stan/math/prim/prob/gamma_lcdf.hpp

@@ -55,7 +55,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
-  // The gradients are technically ill-defined, but treated as zero
+  // The gradients are technically ill-defined
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (y_vec.val(i) == 0) {
       return ops_partials.build(negative_infinity());
@@ -70,7 +70,6 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
     }
 
     const T_partials_return y_dbl = y_vec.val(n);
-    const T_partials_return log_y_dbl = log(y_dbl);
     const T_partials_return alpha_dbl = alpha_vec.val(n);
     const T_partials_return beta_dbl = beta_vec.val(n);
     const T_partials_return log_beta_dbl = log(beta_dbl);
@@ -82,6 +81,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
     P += log_Pn;
 
     if constexpr (is_any_autodiff_v<T_y, T_inv_scale>) {
+      const T_partials_return log_y_dbl = log(y_dbl);
       const T_partials_return d_num
           = (-beta_y_dbl) + (alpha_dbl - 1) * (log_beta_dbl + log_y_dbl);
       const T_partials_return d_den = lgamma(alpha_dbl) + log_Pn;

test/unit/math/prim/prob/gamma_lccdf_test.cpp

@@ -0,0 +1,245 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <limits>
+
+TEST(ProbGamma, lccdf_works) {
+  using stan::math::gamma_lccdf;
+
+  double y = 0.8;
+  double alpha = 1.1;
+  double beta = 2.3;
+
+  EXPECT_NO_THROW(gamma_lccdf(y, alpha, beta));
+}
+
+TEST(ProbGamma, lccdf_zero_y) {
+  using stan::math::gamma_lccdf;
+
+  // When y = 0, LCCDF(0) = log(P(Y > 0)) = log(1) = 0
+  // For continuous distribution, P(Y > 0) = 1
+  double alpha = 1.5;
+  double beta = 2.0;
+
+  double result = gamma_lccdf(0.0, alpha, beta);
+  EXPECT_EQ(result, 0.0);
+}
+
+TEST(ProbGamma, lccdf_large_y) {
+  using stan::math::gamma_lccdf;
+
+  // When y is very large, CDF approaches 1, so LCCDF = log(1-1) = log(0) = -inf
+  double alpha = 1.5;
+  double beta = 2.0;
+  double y = 1e10;
+
+  double result = gamma_lccdf(y, alpha, beta);
+
+  // Should be a very large negative number (approaching -infinity)
+  EXPECT_LT(result, -1000.0);
+}
+
+TEST(ProbGamma, lccdf_infinity_y) {
+  using stan::math::gamma_lccdf;
+  using stan::math::negative_infinity;
+
+  // When y = infinity, LCCDF = log(P(Y > ∞)) = log(0) = -∞
+  double alpha = 1.5;
+  double beta = 2.0;
+  double y = std::numeric_limits<double>::infinity();
+
+  double result = gamma_lccdf(y, alpha, beta);
+  EXPECT_EQ(result, negative_infinity());
+}
+
+TEST(ProbGamma, lccdf_small_alpha_small_y) {
+  using stan::math::gamma_lccdf;
+
+  // Small alpha, small y - numerically challenging
+  double y = 0.001;
+  double alpha = 0.1;
+  double beta = 1.0;
+
+  double result = gamma_lccdf(y, alpha, beta);
+
+  // Should be finite and negative
+  EXPECT_TRUE(std::isfinite(result));
+  EXPECT_LT(result, 0.0);
+}
+
+TEST(ProbGamma, lccdf_large_alpha_large_y) {
+  using stan::math::gamma_lccdf;
+
+  // Large alpha, large y
+  double y = 100.0;
+  double alpha = 50.0;
+  double beta = 0.5;
+
+  double result = gamma_lccdf(y, alpha, beta);
+
+  // Should be finite
+  EXPECT_TRUE(std::isfinite(result));
+}
+
+TEST(ProbGamma, lccdf_alpha_one) {
+  using stan::math::gamma_lccdf;
+  using std::exp;
+  using std::log;
+
+  // When alpha = 1, gamma becomes exponential
+  // For exponential with rate beta: LCCDF(y) = log(1 - (1-exp(-beta*y))) =
+  // log(exp(-beta*y)) = -beta*y
+  double y = 2.0;
+  double alpha = 1.0;
+  double beta = 3.0;
+
+  double result = gamma_lccdf(y, alpha, beta);
+  double expected = -beta * y;  // = -6.0
+
+  EXPECT_NEAR(result, expected, 1e-10);
+}
+
+TEST(ProbGamma, lccdf_various_values) {
+  using stan::math::gamma_lccdf;
+
+  // Test a variety of parameter combinations
+  std::vector<std::tuple<double, double, double>> test_cases = {
+      {0.5, 0.5, 1.0},     // Small y, small alpha
+      {1.0, 1.0, 1.0},     // All ones
+      {2.0, 3.0, 0.5},     // Moderate values
+      {10.0, 2.0, 0.1},    // Large y, small beta
+      {0.1, 10.0, 2.0},    // Small y, large alpha
+      {5.0, 5.0, 1.0},     // Equal alpha and y
+      {0.01, 0.5, 10.0},   // Small y, large beta
+      {100.0, 100.0, 1.0}  // Large matched values
+  };
+
+  for (const auto& test_case : test_cases) {
+    double y = std::get<0>(test_case);
+    double alpha = std::get<1>(test_case);
+    double beta = std::get<2>(test_case);
+
+    double result = gamma_lccdf(y, alpha, beta);
+
+    // All results should be finite and <= 0
+    EXPECT_TRUE(std::isfinite(result))
+        << "Failed for y=" << y << ", alpha=" << alpha << ", beta=" << beta;
+    EXPECT_LE(result, 0.0) << "Failed for y=" << y << ", alpha=" << alpha
+                           << ", beta=" << beta;
+  }
+}
+
+TEST(ProbGamma, lccdf_extreme_small_values) {
+  using stan::math::gamma_lccdf;
+
+  // Very small but non-zero values
+  double y = 1e-10;
+  double alpha = 1e-5;
+  double beta = 1.0;
+
+  double result = gamma_lccdf(y, alpha, beta);
+
+  EXPECT_TRUE(std::isfinite(result));
+}
+
+TEST(ProbGamma, lccdf_extreme_large_alpha) {
+  using stan::math::gamma_lccdf;
+
+  // Very large alpha (approaches normal distribution)
+  double y = 1000.0;
+  double alpha = 1000.0;
+  double beta = 1.0;
+
+  double result = gamma_lccdf(y, alpha, beta);
+
+  EXPECT_TRUE(std::isfinite(result));
+}
+
+TEST(ProbGamma, lccdf_monotonic_in_y) {
+  using stan::math::gamma_lccdf;
+
+  // LCCDF should be monotonically decreasing in y
+  double alpha = 2.0;
+  double beta = 1.5;
+
+  double y1 = 1.0;
+  double y2 = 2.0;
+  double y3 = 3.0;
+
+  double lccdf1 = gamma_lccdf(y1, alpha, beta);
+  double lccdf2 = gamma_lccdf(y2, alpha, beta);
+  double lccdf3 = gamma_lccdf(y3, alpha, beta);
+
+  EXPECT_GT(lccdf1, lccdf2);
+  EXPECT_GT(lccdf2, lccdf3);
+}
+
+TEST(ProbGamma, lccdf_consistency_with_cdf) {
+  using stan::math::gamma_cdf;
+  using stan::math::gamma_lccdf;
+  using std::log;
+
+  // Test that lccdf(y) ≈ log(1 - cdf(y))
+  double y = 1.5;
+  double alpha = 2.5;
+  double beta = 1.8;
+
+  double lccdf_val = gamma_lccdf(y, alpha, beta);
+  double cdf_val = gamma_cdf(y, alpha, beta);
+  double expected = log(1.0 - cdf_val);
+
+  EXPECT_NEAR(lccdf_val, expected, 1e-10);
+}
+
+TEST(ProbGamma, lccdf_numerically_challenging) {
+  using stan::math::gamma_lccdf;
+
+  // Test cases that might cause numerical issues
+  std::vector<std::tuple<double, double, double>> challenging_cases = {
+      {1e-8, 1e-6, 1.0},      // Very small y and alpha
+      {1e-6, 100.0, 1e-3},    // Very small y, large alpha, small beta
+      {1000.0, 0.1, 1e-4},    // Large y, small alpha, very small beta
+      {50.0, 50.0, 1.0},      // Matched moderate values
+      {0.001, 0.001, 100.0},  // Small y and alpha, large beta
+      {1e6, 10.0, 1e-6},      // Very large y, moderate alpha, very small beta
+  };
+
+  for (const auto& test_case : challenging_cases) {
+    double y = std::get<0>(test_case);
+    double alpha = std::get<1>(test_case);
+    double beta = std::get<2>(test_case);
+
+    double result = gamma_lccdf(y, alpha, beta);
+
+    // Should not be NaN
+    EXPECT_FALSE(std::isnan(result))
+        << "NaN for y=" << y << ", alpha=" << alpha << ", beta=" << beta;
+
+    // Should be <= 0 (log of probability)
+    EXPECT_LE(result, 0.0) << "Positive value for y=" << y
+                           << ", alpha=" << alpha << ", beta=" << beta;
+  }
+}
+
+TEST(ProbGamma, lccdf_shape_zero_throws) {
+  using stan::math::gamma_lccdf;
+
+  // alpha (shape) must be positive
+  EXPECT_THROW(gamma_lccdf(1.0, 0.0, 1.0), std::domain_error);
+  EXPECT_THROW(gamma_lccdf(1.0, -1.0, 1.0), std::domain_error);
+}
+
+TEST(ProbGamma, lccdf_rate_zero_throws) {
+  using stan::math::gamma_lccdf;
+
+  // beta (rate) must be positive
+  EXPECT_THROW(gamma_lccdf(1.0, 1.0, 0.0), std::domain_error);
+  EXPECT_THROW(gamma_lccdf(1.0, 1.0, -1.0), std::domain_error);
+}
+
+TEST(ProbGamma, lccdf_negative_y_throws) {
+  using stan::math::gamma_lccdf;
+
+  // y must be non-negative
+  EXPECT_THROW(gamma_lccdf(-1.0, 1.0, 1.0), std::domain_error);
+  EXPECT_THROW(gamma_lccdf(-0.001, 1.0, 1.0), std::domain_error);
+}