1818#include < stan/math/prim/fun/size.hpp>
1919#include < stan/math/prim/fun/size_zero.hpp>
2020#include < stan/math/prim/fun/tgamma.hpp>
21- #include < stan/math/prim/fun/value_of .hpp>
21+ #include < stan/math/prim/fun/value_of_rec .hpp>
2222#include < stan/math/prim/fun/log_gamma_q_dgamma.hpp>
2323#include < stan/math/prim/functor/partials_propagator.hpp>
2424#include < cmath>
@@ -27,12 +27,12 @@ namespace stan {
2727namespace math {
2828
2929template <typename T_y, typename T_shape, typename T_inv_scale>
30- return_type_t <T_y, T_shape, T_inv_scale> gamma_lccdf (const T_y& y,
31- const T_shape& alpha,
32- const T_inv_scale& beta) {
33- using T_partials_return = partials_return_t <T_y, T_shape, T_inv_scale>;
30+ inline return_type_t <T_y, T_shape, T_inv_scale> gamma_lccdf (const T_y& y,
31+ const T_shape& alpha,
32+ const T_inv_scale& beta) {
3433 using std::exp;
3534 using std::log;
35+ using T_partials_return = partials_return_t <T_y, T_shape, T_inv_scale>;
3636 using T_y_ref = ref_type_t <T_y>;
3737 using T_alpha_ref = ref_type_t <T_shape>;
3838 using T_beta_ref = ref_type_t <T_inv_scale>;
@@ -58,7 +58,6 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
5858 scalar_seq_view<T_beta_ref> beta_vec (beta_ref);
5959 const size_t N = max_size (y, alpha, beta);
6060
61- constexpr bool need_y_beta_deriv = !is_constant_all<T_y, T_inv_scale>::value;
6261 constexpr bool any_fvar = is_fvar<scalar_type_t <T_y>>::value
6362 || is_fvar<scalar_type_t <T_shape>>::value
6463 || is_fvar<scalar_type_t <T_inv_scale>>::value;
@@ -91,8 +90,8 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
9190 if constexpr (!any_fvar && is_autodiff_v<T_shape>) {
9291 // var-only path: use log_gamma_q_dgamma which computes both log_q
9392 // and its gradient analytically with double inputs
94- const double beta_y_dbl = value_of ( value_of ( beta_y) );
95- const double alpha_dbl_val = value_of ( value_of ( alpha_dbl) );
93+ const double beta_y_dbl = value_of_rec ( beta_y);
94+ const double alpha_dbl_val = value_of_rec ( alpha_dbl);
9695
9796 if (use_cf) {
9897 auto log_q_result = log_gamma_q_dgamma (alpha_dbl_val, beta_y_dbl);
@@ -105,7 +104,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
105104
106105 // Check if we need to fallback to continued fraction
107106 bool need_cf_fallback
108- = !std::isfinite (value_of ( value_of ( log_Qn) )) || Qn <= 0.0 ;
107+ = !std::isfinite (value_of_rec ( log_Qn)) || Qn <= 0.0 ;
109108 if (need_cf_fallback && beta_y > 0.0 ) {
110109 auto log_q_result = log_gamma_q_dgamma (alpha_dbl_val, beta_y_dbl);
111110 log_Qn = log_q_result.log_q ;
@@ -116,26 +115,26 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
116115 }
117116 } else if constexpr (partials_fvar && is_autodiff_v<T_shape>) {
118117 // fvar path: use unit derivative trick to compute gradients
119- auto alpha_unit = alpha_dbl;
118+ T_partials_return alpha_unit = alpha_dbl;
120119 alpha_unit.d_ = 1 ;
121- auto beta_unit = beta_y;
120+ T_partials_return beta_unit = beta_y;
122121 beta_unit.d_ = 0 ;
123122
124123 if (use_cf) {
125124 log_Qn = internal::log_q_gamma_cf (alpha_dbl, beta_y);
126- auto log_Qn_fvar = internal::log_q_gamma_cf (alpha_unit, beta_unit);
125+ T_partials_return log_Qn_fvar = internal::log_q_gamma_cf (alpha_unit, beta_unit);
127126 dlogQ_dalpha = log_Qn_fvar.d_ ;
128127 } else {
129128 const T_partials_return Pn = gamma_p (alpha_dbl, beta_y);
130129 log_Qn = log1m (Pn);
131130
132- if (!std::isfinite (value_of ( value_of ( log_Qn) )) && beta_y > 0.0 ) {
131+ if (!std::isfinite (value_of_rec ( log_Qn)) && beta_y > 0.0 ) {
133132 // Fallback to continued fraction
134133 log_Qn = internal::log_q_gamma_cf (alpha_dbl, beta_y);
135- auto log_Qn_fvar = internal::log_q_gamma_cf (alpha_unit, beta_unit);
134+ T_partials_return log_Qn_fvar = internal::log_q_gamma_cf (alpha_unit, beta_unit);
136135 dlogQ_dalpha = log_Qn_fvar.d_ ;
137136 } else {
138- auto log_Qn_fvar = log1m (gamma_p (alpha_unit, beta_unit));
137+ T_partials_return log_Qn_fvar = log1m (gamma_p (alpha_unit, beta_unit));
139138 dlogQ_dalpha = log_Qn_fvar.d_ ;
140139 }
141140 }
@@ -147,24 +146,23 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
147146 const T_partials_return Pn = gamma_p (alpha_dbl, beta_y);
148147 log_Qn = log1m (Pn);
149148
150- if (!std::isfinite (value_of ( value_of ( log_Qn) )) && beta_y > 0.0 ) {
149+ if (!std::isfinite (value_of_rec ( log_Qn)) && beta_y > 0.0 ) {
151150 log_Qn = internal::log_q_gamma_cf (alpha_dbl, beta_y);
152151 }
153152 }
154153 }
155- if (!std::isfinite (value_of ( value_of ( log_Qn) ))) {
154+ if (!std::isfinite (value_of_rec ( log_Qn))) {
156155 return ops_partials.build (negative_infinity ());
157156 }
158157 P += log_Qn;
159158
160- if constexpr (need_y_beta_deriv) {
159+ constexpr bool need_y_beta_deriv = !is_constant_all<T_y, T_inv_scale>::value;
160+ if constexpr (is_autodiff_v<T_y> || is_autodiff_v<T_inv_scale>) {
161161 const T_partials_return log_y = log (y_dbl);
162- const T_partials_return log_beta = log (beta_dbl);
163- const T_partials_return lgamma_alpha = lgamma (alpha_dbl);
164162 const T_partials_return alpha_minus_one = fma (alpha_dbl, log_y, -log_y);
165163
166164 const T_partials_return log_pdf
167- = alpha_dbl * log_beta - lgamma_alpha + alpha_minus_one - beta_y;
165+ = alpha_dbl * log (beta_dbl) - lgamma (alpha_dbl) + alpha_minus_one - beta_y;
168166
169167 const T_partials_return hazard = exp (log_pdf - log_Qn); // f/Q
170168
0 commit comments