2626namespace stan {
2727namespace math {
2828namespace internal {
29- template <typename T>
30- struct Q_eval {
31- T log_Q{0.0 };
32- T dlogQ_dalpha{0.0 };
33- bool ok{false };
34- };
35-
36- /* *
37- * Computes log q and d(log q) / d(alpha) using continued fraction.
38- */
39- template <typename T, typename T_shape,
40- bool any_fvar, bool partials_fvar>
41- static inline Q_eval<T> eval_q_cf (const T& alpha,
42- const T& beta_y) {
43- Q_eval<T> out;
44- if constexpr (!any_fvar && is_autodiff_v<T_shape>) {
45- auto log_q_result = log_gamma_q_dgamma (value_of_rec (alpha), value_of_rec (beta_y));
46- out.log_Q = log_q_result.log_q ;
47- out.dlogQ_dalpha = log_q_result.dlog_q_da ;
48- } else {
49- out.log_Q = internal::log_q_gamma_cf (alpha, beta_y);
50- if constexpr (is_autodiff_v<T_shape>) {
51- if constexpr (!partials_fvar) {
52- out.dlogQ_dalpha
53- = grad_reg_inc_gamma (alpha, beta_y, tgamma (alpha),
54- digamma (alpha)) / exp (out.log_Q );
55- } else {
56- T alpha_unit = alpha;
57- alpha_unit.d_ = 1 ;
58- T beta_y_unit = beta_y;
59- beta_y_unit.d_ = 0 ;
60- T log_Q_fvar = internal::log_q_gamma_cf (alpha_unit, beta_y_unit);
61- out.dlogQ_dalpha = log_Q_fvar.d_ ;
62- }
63- }
64- }
65-
66- out.ok = std::isfinite (value_of_rec (out.log_Q ));
67- return out;
68- }
69-
70- /* *
71- * Computes log q and d(log q) / d(alpha) using log1m.
72- */
73- template <typename T, typename T_shape,
74- bool partials_fvar>
75- static inline Q_eval<T> eval_q_log1m (const T& alpha,
76- const T& beta_y) {
77- Q_eval<T> out;
78- out.log_Q = log1m (gamma_p (alpha, beta_y));
79-
80- if (!std::isfinite (value_of_rec (out.log_Q ))) {
81- out.ok = false ;
82- return out;
83- }
84-
29+ template <typename T>
30+ struct Q_eval {
31+ T log_Q{0.0 };
32+ T dlogQ_dalpha{0.0 };
33+ bool ok{false };
34+ };
35+
36+ /* *
37+ * Computes log q and d(log q) / d(alpha) using continued fraction.
38+ */
39+ template <typename T, typename T_shape, bool any_fvar, bool partials_fvar>
40+ static inline Q_eval<T> eval_q_cf (const T& alpha, const T& beta_y) {
41+ Q_eval<T> out;
42+ if constexpr (!any_fvar && is_autodiff_v<T_shape>) {
43+ auto log_q_result
44+ = log_gamma_q_dgamma (value_of_rec (alpha), value_of_rec (beta_y));
45+ out.log_Q = log_q_result.log_q ;
46+ out.dlogQ_dalpha = log_q_result.dlog_q_da ;
47+ } else {
48+ out.log_Q = internal::log_q_gamma_cf (alpha, beta_y);
8549 if constexpr (is_autodiff_v<T_shape>) {
86- if constexpr (partials_fvar) {
50+ if constexpr (!partials_fvar) {
51+ out.dlogQ_dalpha
52+ = grad_reg_inc_gamma (alpha, beta_y, tgamma (alpha), digamma (alpha))
53+ / exp (out.log_Q );
54+ } else {
8755 T alpha_unit = alpha;
8856 alpha_unit.d_ = 1 ;
89- T beta_unit = beta_y;
90- beta_unit .d_ = 0 ;
91- T log_Q_fvar = log1m ( gamma_p ( alpha_unit, beta_unit) );
57+ T beta_y_unit = beta_y;
58+ beta_y_unit .d_ = 0 ;
59+ T log_Q_fvar = internal::log_q_gamma_cf ( alpha_unit, beta_y_unit );
9260 out.dlogQ_dalpha = log_Q_fvar.d_ ;
93- } else {
94- out.dlogQ_dalpha = -grad_reg_lower_inc_gamma (alpha, beta_y) / exp (out.log_Q );
9561 }
9662 }
63+ }
64+
65+ out.ok = std::isfinite (value_of_rec (out.log_Q ));
66+ return out;
67+ }
9768
98- out.ok = true ;
69+ /* *
70+ * Computes log q and d(log q) / d(alpha) using log1m.
71+ */
72+ template <typename T, typename T_shape, bool partials_fvar>
73+ static inline Q_eval<T> eval_q_log1m (const T& alpha, const T& beta_y) {
74+ Q_eval<T> out;
75+ out.log_Q = log1m (gamma_p (alpha, beta_y));
76+
77+ if (!std::isfinite (value_of_rec (out.log_Q ))) {
78+ out.ok = false ;
9979 return out;
10080 }
81+
82+ if constexpr (is_autodiff_v<T_shape>) {
83+ if constexpr (partials_fvar) {
84+ T alpha_unit = alpha;
85+ alpha_unit.d_ = 1 ;
86+ T beta_unit = beta_y;
87+ beta_unit.d_ = 0 ;
88+ T log_Q_fvar = log1m (gamma_p (alpha_unit, beta_unit));
89+ out.dlogQ_dalpha = log_Q_fvar.d_ ;
90+ } else {
91+ out.dlogQ_dalpha
92+ = -grad_reg_lower_inc_gamma (alpha, beta_y) / exp (out.log_Q );
93+ }
94+ }
95+
96+ out.ok = true ;
97+ return out;
10198}
99+ } // namespace internal
102100
103101template <typename T_y, typename T_shape, typename T_inv_scale>
104- inline return_type_t <T_y, T_shape, T_inv_scale> gamma_lccdf (const T_y& y,
105- const T_shape& alpha,
106- const T_inv_scale& beta) {
102+ inline return_type_t <T_y, T_shape, T_inv_scale> gamma_lccdf (
103+ const T_y& y, const T_shape& alpha, const T_inv_scale& beta) {
107104 using std::exp;
108105 using std::log;
109106 using T_partials_return = partials_return_t <T_y, T_shape, T_inv_scale>;
@@ -159,16 +156,17 @@ inline return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
159156 const bool use_continued_fraction = beta_y > alpha_dbl + 1.0 ;
160157 internal::Q_eval<T_partials_return> result;
161158 if (use_continued_fraction) {
162- result = internal::eval_q_cf<T_partials_return, T_shape,
163- any_fvar, partials_fvar>(alpha_dbl, beta_y);
164- } else {
165- result = internal::eval_q_log1m<T_partials_return, T_shape,
159+ result = internal::eval_q_cf<T_partials_return, T_shape, any_fvar,
166160 partials_fvar>(alpha_dbl, beta_y);
161+ } else {
162+ result
163+ = internal::eval_q_log1m<T_partials_return, T_shape, partials_fvar>(
164+ alpha_dbl, beta_y);
167165
168166 if (!result.ok && beta_y > 0.0 ) {
169167 // Fallback to continued fraction if log1m fails
170- result = internal::eval_q_cf<T_partials_return, T_shape,
171- any_fvar, partials_fvar>(alpha_dbl, beta_y);
168+ result = internal::eval_q_cf<T_partials_return, T_shape, any_fvar,
169+ partials_fvar>(alpha_dbl, beta_y);
172170 }
173171 }
174172 if (!result.ok ) {
@@ -181,8 +179,9 @@ inline return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
181179 const T_partials_return log_y = log (y_dbl);
182180 const T_partials_return alpha_minus_one = fma (alpha_dbl, log_y, -log_y);
183181
184- const T_partials_return log_pdf
185- = alpha_dbl * log (beta_dbl) - lgamma (alpha_dbl) + alpha_minus_one - beta_y;
182+ const T_partials_return log_pdf = alpha_dbl * log (beta_dbl)
183+ - lgamma (alpha_dbl) + alpha_minus_one
184+ - beta_y;
186185
187186 const T_partials_return hazard = exp (log_pdf - result.log_Q ); // f/Q
188187
0 commit comments