Automatic Differentiation
 
Loading...
Searching...
No Matches
beta_cdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_BETA_CDF_HPP
2#define STAN_MATH_PRIM_PROB_BETA_CDF_HPP
3
18#include <cmath>
19
20namespace stan {
21namespace math {
22
35template <typename T_y, typename T_scale_succ, typename T_scale_fail>
37 const T_y& y, const T_scale_succ& alpha, const T_scale_fail& beta) {
39 using T_y_ref = ref_type_t<T_y>;
40 using T_alpha_ref = ref_type_t<T_scale_succ>;
41 using T_beta_ref = ref_type_t<T_scale_fail>;
42 static constexpr const char* function = "beta_cdf";
43 check_consistent_sizes(function, "Random variable", y,
44 "First shape parameter", alpha,
45 "Second shape parameter", beta);
46 if (size_zero(y, alpha, beta)) {
47 return 1.0;
48 }
49
50 T_y_ref y_ref = y;
51 T_alpha_ref alpha_ref = alpha;
52 T_beta_ref beta_ref = beta;
53 check_positive_finite(function, "First shape parameter", alpha_ref);
54 check_positive_finite(function, "Second shape parameter", beta_ref);
55 check_bounded(function, "Random variable", value_of(y_ref), 0, 1);
56
57 T_partials_return P(1.0);
58 auto ops_partials = make_partials_propagator(y_ref, alpha_ref, beta_ref);
59 scalar_seq_view<T_y_ref> y_vec(y_ref);
60 scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
61 scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
62 size_t size_alpha = stan::math::size(alpha);
63 size_t size_beta = stan::math::size(beta);
64 size_t size_alpha_beta = max_size(alpha, beta);
65 size_t N = max_size(y, alpha, beta);
66
67 // Explicit return for extreme values
68 // The gradients are technically ill-defined, but treated as zero
69 for (size_t i = 0; i < stan::math::size(y); i++) {
70 if (y_vec.val(i) <= 0) {
71 return ops_partials.build(0.0);
72 }
73 }
74
76 T_scale_succ>
77 digamma_alpha(size_alpha);
79 for (size_t n = 0; n < size_alpha; n++) {
80 digamma_alpha[n] = digamma(alpha_vec.val(n));
81 }
82 }
83
85 T_scale_fail>
86 digamma_beta(size_beta);
88 for (size_t n = 0; n < size_beta; n++) {
89 digamma_beta[n] = digamma(beta_vec.val(n));
90 }
91 }
92
94 T_partials_return, T_scale_succ, T_scale_fail>
95 digamma_sum(size_alpha_beta);
97 for (size_t n = 0; n < size_alpha_beta; n++) {
98 digamma_sum[n] = digamma(alpha_vec.val(n) + beta_vec.val(n));
99 }
100 }
101
102 for (size_t n = 0; n < N; n++) {
103 const T_partials_return y_dbl = y_vec.val(n);
104
105 // Explicit results for extreme values
106 // The gradients are technically ill-defined, but treated as zero
107 if (y_dbl >= 1.0) {
108 continue;
109 }
110
111 const T_partials_return alpha_dbl = alpha_vec.val(n);
112 const T_partials_return beta_dbl = beta_vec.val(n);
113 const T_partials_return Pn = inc_beta(alpha_dbl, beta_dbl, y_dbl);
114 const T_partials_return inv_Pn
116
117 P *= Pn;
118
120 partials<0>(ops_partials)[n]
121 += inc_beta_ddz(alpha_dbl, beta_dbl, y_dbl) * inv_Pn;
122 }
123
125 partials<1>(ops_partials)[n]
126 += inc_beta_dda(alpha_dbl, beta_dbl, y_dbl, digamma_alpha[n],
127 digamma_sum[n])
128 * inv_Pn;
129 }
131 partials<2>(ops_partials)[n]
132 += inc_beta_ddb(alpha_dbl, beta_dbl, y_dbl, digamma_beta[n],
133 digamma_sum[n])
134 * inv_Pn;
135 }
136 }
137
139 for (size_t n = 0; n < stan::math::size(y); ++n) {
140 partials<0>(ops_partials)[n] *= P;
141 }
142 }
144 for (size_t n = 0; n < stan::math::size(alpha); ++n) {
145 partials<1>(ops_partials)[n] *= P;
146 }
147 }
149 for (size_t n = 0; n < stan::math::size(beta); ++n) {
150 partials<2>(ops_partials)[n] *= P;
151 }
152 }
153
154 return ops_partials.build(P);
155}
156
157} // namespace math
158} // namespace stan
159#endif
VectorBuilder allocates type T1 values to be used as intermediate values.
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
return_type_t< T_y, T_scale_succ, T_scale_fail > beta_cdf(const T_y &y, const T_scale_succ &alpha, const T_scale_fail &beta)
Calculates the beta cumulative distribution function for the given variate and scale variables.
Definition beta_cdf.hpp:36
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
T inc_beta_ddz(T a, T b, T z)
Returns the partial derivative of the regularized incomplete beta function, I_{z}(a,...
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Check if the value is between the low and high values, inclusively.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
T inc_beta_dda(T a, T b, T z, T digamma_a, T digamma_ab)
Returns the partial derivative of the regularized incomplete beta function, I_{z}(a,...
fvar< T > inc_beta(const fvar< T > &a, const fvar< T > &b, const fvar< T > &x)
Definition inc_beta.hpp:19
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
fvar< T > inv(const fvar< T > &x)
Definition inv.hpp:12
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.
fvar< T > digamma(const fvar< T > &x)
Return the derivative of the log gamma function at the specified argument.
Definition digamma.hpp:23
T inc_beta_ddb(T a, T b, T z, T digamma_b, T digamma_ab)
Returns the partial derivative of the regularized incomplete beta function, I_{z}(a,...
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...