Automatic Differentiation
 
Loading...
Searching...
No Matches
gamma_lcdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_GAMMA_LCDF_HPP
2#define STAN_MATH_PRIM_PROB_GAMMA_LCDF_HPP
3
19#include <cmath>
20
21namespace stan {
22namespace math {
23
24template <typename T_y, typename T_shape, typename T_inv_scale>
26 const T_shape& alpha,
27 const T_inv_scale& beta) {
28 using T_partials_return = partials_return_t<T_y, T_shape, T_inv_scale>;
29 using std::exp;
30 using std::log;
31 using std::pow;
32 using T_y_ref = ref_type_t<T_y>;
33 using T_alpha_ref = ref_type_t<T_shape>;
34 using T_beta_ref = ref_type_t<T_inv_scale>;
35 static constexpr const char* function = "gamma_lcdf";
36 check_consistent_sizes(function, "Random variable", y, "Shape parameter",
37 alpha, "Inverse scale parameter", beta);
38 T_y_ref y_ref = y;
39 T_alpha_ref alpha_ref = alpha;
40 T_beta_ref beta_ref = beta;
41 check_positive_finite(function, "Shape parameter", alpha_ref);
42 check_positive_finite(function, "Inverse scale parameter", beta_ref);
43 check_nonnegative(function, "Random variable", y_ref);
44
45 if (size_zero(y, alpha, beta)) {
46 return 0;
47 }
48
49 T_partials_return P(0.0);
50 auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref);
51
52 scalar_seq_view<T_y_ref> y_vec(y_ref);
53 scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
54 scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
55 size_t N = max_size(y, alpha, beta);
56
57 // Explicit return for extreme values
58 // The gradients are technically ill-defined, but treated as zero
59 for (size_t i = 0; i < stan::math::size(y); i++) {
60 if (y_vec.val(i) == 0) {
61 return ops_partials.build(negative_infinity());
62 }
63 }
64
65 for (size_t n = 0; n < N; n++) {
66 // Explicit results for extreme values
67 // The gradients are technically ill-defined, but treated as zero
68 if (y_vec.val(n) == INFTY) {
69 return ops_partials.build(0.0);
70 }
71
72 const T_partials_return y_dbl = y_vec.val(n);
73 const T_partials_return log_y_dbl = log(y_dbl);
74 const T_partials_return alpha_dbl = alpha_vec.val(n);
75 const T_partials_return beta_dbl = beta_vec.val(n);
76 const T_partials_return log_beta_dbl = log(beta_dbl);
77 const T_partials_return beta_y_dbl = beta_dbl * y_dbl;
78
79 const T_partials_return Pn = gamma_p(alpha_dbl, beta_y_dbl);
80 const T_partials_return log_Pn = log(Pn);
81
82 P += log_Pn;
83
85 const T_partials_return d_num
86 = (-beta_y_dbl) + (alpha_dbl - 1) * (log_beta_dbl + log_y_dbl);
87 const T_partials_return d_den = lgamma(alpha_dbl) + log_Pn;
88 const T_partials_return d = exp(d_num - d_den);
89
91 partials<0>(ops_partials)[n] += beta_dbl * d;
92 }
94 partials<2>(ops_partials)[n] += y_dbl * d;
95 }
96 }
98 partials<1>(ops_partials)[n]
99 += grad_reg_lower_inc_gamma(alpha_dbl, beta_y_dbl) / Pn;
100 }
101 }
102 return ops_partials.build(P);
103}
104
105} // namespace math
106} // namespace stan
107#endif
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
static constexpr double negative_infinity()
Return negative infinity.
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
return_type_t< T_y, T_shape, T_inv_scale > gamma_lcdf(const T_y &y, const T_shape &alpha, const T_inv_scale &beta)
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
fvar< T > lgamma(const fvar< T > &x)
Return the natural logarithm of the gamma function applied to the specified argument.
Definition lgamma.hpp:21
fvar< T > gamma_p(const fvar< T > &x1, const fvar< T > &x2)
Definition gamma_p.hpp:16
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
fvar< T > beta(const fvar< T > &x1, const fvar< T > &x2)
Return fvar with the beta function applied to the specified arguments and its gradient.
Definition beta.hpp:51
return_type_t< T1, T2 > grad_reg_lower_inc_gamma(const T1 &a, const T2 &z, double precision=1e-10, int max_steps=1e5)
Computes the gradient of the lower regularized incomplete gamma function.
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
typename partials_return_type< Args... >::type partials_return_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Extends std::true_type when instantiated with zero or more template parameters, all of which extend t...