Automatic Differentiation
 
Loading...
Searching...
No Matches
beta_binomial_lpmf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_BETA_BINOMIAL_LPMF_HPP
2#define STAN_MATH_PRIM_PROB_BETA_BINOMIAL_LPMF_HPP
3
18
19namespace stan {
20namespace math {
21
39template <bool propto, typename T_n, typename T_N, typename T_size1,
40 typename T_size2,
42 T_n, T_N, T_size1, T_size2>* = nullptr>
44 const T_size1& alpha,
45 const T_size2& beta) {
46 using T_partials_return = partials_return_t<T_size1, T_size2>;
47 using T_N_ref = ref_type_t<T_N>;
48 using T_alpha_ref = ref_type_t<T_size1>;
49 using T_beta_ref = ref_type_t<T_size2>;
50 static constexpr const char* function = "beta_binomial_lpmf";
51 check_consistent_sizes(function, "Successes variable", n,
52 "Population size parameter", N,
53 "First prior sample size parameter", alpha,
54 "Second prior sample size parameter", beta);
55 if (size_zero(n, N, alpha, beta)) {
56 return 0.0;
57 }
58
59 T_N_ref N_ref = N;
60 T_alpha_ref alpha_ref = alpha;
61 T_beta_ref beta_ref = beta;
62 check_nonnegative(function, "Population size parameter", N_ref);
63 check_positive_finite(function, "First prior sample size parameter",
64 alpha_ref);
65 check_positive_finite(function, "Second prior sample size parameter",
66 beta_ref);
67
69 return 0.0;
70 }
71
72 T_partials_return logp(0.0);
73 auto ops_partials = make_partials_propagator(alpha_ref, beta_ref);
74
75 scalar_seq_view<T_n> n_vec(n);
76 scalar_seq_view<T_N_ref> N_vec(N_ref);
77 scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
78 scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
79 size_t size_alpha = stan::math::size(alpha);
80 size_t size_beta = stan::math::size(beta);
81 size_t size_n_N = max_size(n, N);
82 size_t size_alpha_beta = max_size(alpha, beta);
83 size_t max_size_seq_view = max_size(n, N, alpha, beta);
84
85 for (size_t i = 0; i < max_size_seq_view; i++) {
86 if (n_vec[i] < 0 || n_vec[i] > N_vec[i]) {
87 return ops_partials.build(LOG_ZERO);
88 }
89 }
90
91 VectorBuilder<include_summand<propto>::value, T_partials_return, T_n, T_N>
92 normalizing_constant(size_n_N);
93 for (size_t i = 0; i < size_n_N; i++)
95 normalizing_constant[i] = binomial_coefficient_log(N_vec[i], n_vec[i]);
96
98 size_alpha_beta);
99 for (size_t i = 0; i < size_alpha_beta; i++) {
100 lbeta_denominator[i] = lbeta(alpha_vec.val(i), beta_vec.val(i));
101 }
102
104 max_size_seq_view);
105 for (size_t i = 0; i < max_size_seq_view; i++) {
106 lbeta_diff[i] = lbeta(n_vec[i] + alpha_vec.val(i),
107 N_vec[i] - n_vec[i] + beta_vec.val(i))
108 - lbeta_denominator[i];
109 }
110
112 T_size1>
113 digamma_n_plus_alpha(max_size(n, alpha));
115 for (size_t i = 0; i < max_size(n, alpha); i++) {
116 digamma_n_plus_alpha[i] = digamma(n_vec.val(i) + alpha_vec.val(i));
117 }
118 }
119
121 T_size1, T_size2>
122 digamma_alpha_plus_beta(size_alpha_beta);
124 for (size_t i = 0; i < size_alpha_beta; i++) {
125 digamma_alpha_plus_beta[i] = digamma(alpha_vec.val(i) + beta_vec.val(i));
126 }
127 }
128
130 T_N, T_size1, T_size2>
131 digamma_diff(max_size(N, alpha, beta));
133 for (size_t i = 0; i < max_size(N, alpha, beta); i++) {
134 digamma_diff[i]
135 = digamma_alpha_plus_beta[i]
136 - digamma(N_vec.val(i) + alpha_vec.val(i) + beta_vec.val(i));
137 }
138 }
139
140 VectorBuilder<!is_constant_all<T_size1>::value, T_partials_return, T_size1>
141 digamma_alpha(size_alpha);
142 for (size_t i = 0; i < size_alpha; i++)
144 digamma_alpha[i] = digamma(alpha_vec.val(i));
145
146 VectorBuilder<!is_constant_all<T_size2>::value, T_partials_return, T_size2>
147 digamma_beta(size_beta);
148 for (size_t i = 0; i < size_beta; i++)
150 digamma_beta[i] = digamma(beta_vec.val(i));
151
152 for (size_t i = 0; i < max_size_seq_view; i++) {
154 logp += normalizing_constant[i];
155 logp += lbeta_diff[i];
156
158 partials<0>(ops_partials)[i]
159 += digamma_n_plus_alpha[i] + digamma_diff[i] - digamma_alpha[i];
161 partials<1>(ops_partials)[i]
162 += digamma(N_vec.val(i) - n_vec.val(i) + beta_vec.val(i))
163 + digamma_diff[i] - digamma_beta[i];
164 }
165 return ops_partials.build(logp);
166}
167
168template <typename T_n, typename T_N, typename T_size1, typename T_size2>
170 const T_size1& alpha,
171 const T_size2& beta) {
172 return beta_binomial_lpmf<false>(n, N, alpha, beta);
173}
174
175} // namespace math
176} // namespace stan
177#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...
require_all_not_t< is_nonscalar_prim_or_rev_kernel_expression< std::decay_t< Types > >... > require_all_not_nonscalar_prim_or_rev_kernel_expression_t
Require none of the types satisfy is_nonscalar_prim_or_rev_kernel_expression.
binomial_coefficient_log_< as_operation_cl_t< T1 >, as_operation_cl_t< T2 > > binomial_coefficient_log(T1 &&a, T2 &&b)
return_type_t< T_n_cl, T_size1_cl, T_size2_cl > beta_binomial_lpmf(const T_n_cl &n, const T_N_cl N, const T_size1_cl &alpha, const T_size2_cl &beta)
Returns the log PMF of the Beta-Binomial distribution with given population size, prior success,...
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 LOG_ZERO
The natural logarithm of 0, .
Definition constants.hpp:68
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
fvar< T > lbeta(const fvar< T > &x1, const fvar< T > &x2)
Definition lbeta.hpp:14
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
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
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
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...
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...