Skip to content

Commit 44ba284

Browse files
spinkneysyclik
spinkney
authored and
syclik
committed
remove internal function in lccdf and add two tests for expanded range
1 parent bb7391c commit 44ba284

4 files changed

Lines changed: 107 additions & 94 deletions

File tree

‎stan/math/prim/fun/log_gamma_q_dgamma.hpp‎

Lines changed: 62 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,66 @@ struct log_gamma_q_result {
2929
T dlog_q_da; ///< d/da log(Q(a,z))
3030
};
3131

32+
namespace internal {
33+
34+
/**
35+
* Compute log(Q(a,z)) using continued fraction expansion for upper incomplete
36+
* gamma function.
37+
*
38+
* @tparam T_a Type of shape parameter a (double or fvar types)
39+
* @tparam T_z Type of value parameter z (double or fvar types)
40+
* @param a Shape parameter
41+
* @param z Value at which to evaluate
42+
* @param max_steps Maximum number of continued fraction iterations
43+
* @param precision Convergence threshold
44+
* @return log(Q(a,z)) with same type as T_a and T_z
45+
*/
46+
template <typename T_a, typename T_z>
47+
inline auto log_q_gamma_cf(const T_a& a, const T_z& z, int max_steps = 250,
48+
double precision = 1e-16) {
49+
using stan::math::lgamma;
50+
using stan::math::log;
51+
using stan::math::value_of;
52+
using std::fabs;
53+
using T_return = return_type_t<T_a, T_z>;
54+
55+
const T_return a_ret = a;
56+
const T_return z_ret = z;
57+
const auto log_prefactor = a_ret * log(z_ret) - z_ret - lgamma(a_ret);
58+
59+
auto b = z_ret + 1.0 - a_ret;
60+
auto C = (fabs(value_of(b)) >= EPSILON) ? b : T_return(EPSILON);
61+
auto D = T_return(0.0);
62+
auto f = C;
63+
64+
for (int i = 1; i <= max_steps; ++i) {
65+
auto an = -i * (i - a_ret);
66+
b += 2.0;
67+
68+
D = b + an * D;
69+
if (fabs(value_of(D)) < EPSILON) {
70+
D = T_return(EPSILON);
71+
}
72+
C = b + an / C;
73+
if (fabs(value_of(C)) < EPSILON) {
74+
C = T_return(EPSILON);
75+
}
76+
77+
D = 1.0 / D;
78+
auto delta = C * D;
79+
f *= delta;
80+
81+
const double delta_m1 = value_of(fabs(value_of(delta) - 1.0));
82+
if (delta_m1 < precision) {
83+
break;
84+
}
85+
}
86+
87+
return log_prefactor - log(f);
88+
}
89+
90+
} // namespace internal
91+
3292
/**
3393
* Compute log(Q(a,z)) and its gradient with respect to a using continued
3494
* fraction expansion, where Q(a,z) = Gamma(a,z) / Gamma(a) is the regularized
@@ -50,7 +110,6 @@ template <typename T_a, typename T_z>
50110
inline log_gamma_q_result<return_type_t<T_a, T_z>> log_gamma_q_dgamma(
51111
const T_a& a, const T_z& z, int max_steps = 250, double precision = 1e-16) {
52112
using std::exp;
53-
using std::fabs;
54113
using std::log;
55114
using T_return = return_type_t<T_a, T_z>;
56115

@@ -61,39 +120,8 @@ inline log_gamma_q_result<return_type_t<T_a, T_z>> log_gamma_q_dgamma(
61120

62121
// For z > a + 1, use continued fraction for better numerical stability
63122
if (z_dbl > a_dbl + 1.0) {
64-
// Continued fraction for Q(a,z) in log space
65-
// log(Q(a,z)) = log_prefactor - log(continued_fraction)
66-
const double log_prefactor = a_dbl * log(z_dbl) - z_dbl - lgamma(a_dbl);
67-
68-
double b = z_dbl + 1.0 - a_dbl;
69-
double C = (fabs(b) >= EPSILON) ? b : EPSILON;
70-
double D = 0.0;
71-
double f = C;
72-
73-
for (int i = 1; i <= max_steps; ++i) {
74-
const double an = -i * (i - a_dbl);
75-
b += 2.0;
76-
77-
D = b + an * D;
78-
if (fabs(D) < EPSILON) {
79-
D = EPSILON;
80-
}
81-
C = b + an / C;
82-
if (fabs(C) < EPSILON) {
83-
C = EPSILON;
84-
}
85-
86-
D = 1.0 / D;
87-
const double delta = C * D;
88-
f *= delta;
89-
90-
const double delta_m1 = fabs(delta - 1.0);
91-
if (delta_m1 < precision) {
92-
break;
93-
}
94-
}
95-
96-
result.log_q = log_prefactor - log(f);
123+
result.log_q
124+
= internal::log_q_gamma_cf(a_dbl, z_dbl, max_steps, precision);
97125

98126
// For gradient, use: d/da log(Q) = (1/Q) * dQ/da
99127
// grad_reg_inc_gamma computes dQ/da

‎stan/math/prim/prob/gamma_lccdf.hpp‎

Lines changed: 0 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -26,66 +26,6 @@
2626
namespace stan {
2727
namespace math {
2828

29-
namespace internal {
30-
31-
/**
32-
* Compute log(Q(a,x)) using continued fraction expansion for upper incomplete
33-
* gamma function. When used with fvar types, automatically computes
34-
* derivatives.
35-
*
36-
* @tparam T_a Type of shape parameter a (double or fvar types)
37-
* @param a Shape parameter
38-
* @param x Value at which to evaluate
39-
* @param max_steps Maximum number of continued fraction iterations
40-
* @param precision Convergence threshold
41-
* @return log(Q(a,x)) with same type as T_a
42-
*/
43-
template <typename T_a, typename T_x>
44-
inline auto log_q_gamma_cf(const T_a& a, const T_x& x, int max_steps = 250,
45-
double precision = 1e-16) {
46-
using stan::math::lgamma;
47-
using stan::math::log;
48-
using stan::math::value_of;
49-
using std::fabs;
50-
using T_return = return_type_t<T_a, T_x>;
51-
52-
const T_return a_ret = a;
53-
const T_return x_ret = x;
54-
const auto log_prefactor = a_ret * log(x_ret) - x_ret - lgamma(a_ret);
55-
56-
auto b = x_ret + 1.0 - a_ret;
57-
auto C = (fabs(value_of(b)) >= EPSILON) ? b : T_return(EPSILON);
58-
auto D = T_return(0.0);
59-
auto f = C;
60-
61-
for (int i = 1; i <= max_steps; ++i) {
62-
auto an = -i * (i - a_ret);
63-
b += 2.0;
64-
65-
D = b + an * D;
66-
if (fabs(value_of(D)) < EPSILON) {
67-
D = T_a(EPSILON);
68-
}
69-
C = b + an / C;
70-
if (fabs(value_of(C)) < EPSILON) {
71-
C = T_a(EPSILON);
72-
}
73-
74-
D = 1.0 / D;
75-
auto delta = C * D;
76-
f *= delta;
77-
78-
const double delta_m1 = value_of(fabs(value_of(delta) - 1.0));
79-
if (delta_m1 < precision) {
80-
break;
81-
}
82-
}
83-
84-
return log_prefactor - log(f);
85-
}
86-
87-
} // namespace internal
88-
8929
template <typename T_y, typename T_shape, typename T_inv_scale>
9030
return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
9131
const T_shape& alpha,

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

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,26 @@ TEST(ProbGamma, lccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) {
9191
EXPECT_DOUBLE_EQ(new_val, expected);
9292
}
9393

94+
TEST(ProbGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) {
95+
using stan::math::gamma_lcdf;
96+
using stan::math::gamma_lccdf;
97+
using stan::math::log1m_exp;
98+
using stan::math::negative_infinity;
99+
100+
double y = 20000.0;
101+
double alpha = 800.0;
102+
double beta = 0.1;
103+
104+
double lcdf = gamma_lcdf(y, alpha, beta);
105+
double log1m_lcdf = log1m_exp(lcdf);
106+
double lccdf = gamma_lccdf(y, alpha, beta);
107+
108+
EXPECT_EQ(lcdf, 0.0);
109+
EXPECT_EQ(log1m_lcdf, negative_infinity());
110+
EXPECT_TRUE(std::isfinite(lccdf));
111+
EXPECT_LT(lccdf, 0.0);
112+
}
113+
94114
TEST(ProbGamma, lccdf_large_alpha_large_y) {
95115
using stan::math::gamma_lccdf;
96116

‎test/unit/math/rev/prob/gamma_lccdf_test.cpp‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -272,6 +272,31 @@ TEST(ProbDistributionsGamma,
272272
EXPECT_LE(grads[0], 0.0);
273273
}
274274

275+
TEST(ProbDistributionsGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) {
276+
using stan::math::gamma_lcdf;
277+
using stan::math::gamma_lccdf;
278+
using stan::math::log1m_exp;
279+
using stan::math::negative_infinity;
280+
using stan::math::var;
281+
282+
double y_d = 20000.0;
283+
double alpha_d = 800.0;
284+
double beta_d = 0.1;
285+
286+
double lcdf = gamma_lcdf(y_d, alpha_d, beta_d);
287+
double log1m_lcdf = log1m_exp(lcdf);
288+
289+
var y_v = y_d;
290+
var alpha_v = alpha_d;
291+
var beta_v = beta_d;
292+
var lccdf_var = gamma_lccdf(y_v, alpha_v, beta_v);
293+
294+
EXPECT_EQ(lcdf, 0.0);
295+
EXPECT_EQ(log1m_lcdf, negative_infinity());
296+
EXPECT_TRUE(std::isfinite(lccdf_var.val()));
297+
EXPECT_LT(lccdf_var.val(), 0.0);
298+
}
299+
275300
TEST(ProbDistributionsGamma, lccdf_extreme_values_large) {
276301
using stan::math::gamma_lccdf;
277302
using stan::math::var;

0 commit comments

Comments
 (0)