Automatic Differentiation
 
Loading...
Searching...
No Matches
binomial_lpmf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_BINOMIAL_LPMF_HPP
2#define STAN_MATH_PRIM_PROB_BINOMIAL_LPMF_HPP
3
15
16namespace stan {
17namespace math {
18
36template <bool propto, typename T_n, typename T_N, typename T_prob,
38 T_n, T_N, T_prob>* = nullptr>
39return_type_t<T_prob> binomial_lpmf(const T_n& n, const T_N& N,
40 const T_prob& theta) {
41 using T_partials_return = partials_return_t<T_n, T_N, T_prob>;
42 using T_n_ref = ref_type_t<T_n>;
43 using T_N_ref = ref_type_t<T_N>;
44 using T_theta_ref = ref_type_t<T_prob>;
45 static constexpr const char* function = "binomial_lpmf";
46 check_consistent_sizes(function, "Successes variable", n,
47 "Population size parameter", N,
48 "Probability parameter", theta);
49
50 T_n_ref n_ref = n;
51 T_N_ref N_ref = N;
52 T_theta_ref theta_ref = theta;
53
54 check_bounded(function, "Successes variable", value_of(n_ref), 0, N_ref);
55 check_nonnegative(function, "Population size parameter", N_ref);
56 check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0,
57 1.0);
58
59 if (size_zero(n, N, theta)) {
60 return 0.0;
61 }
63 return 0.0;
64 }
65
66 T_partials_return logp = 0;
67 auto ops_partials = make_partials_propagator(theta_ref);
68
69 scalar_seq_view<T_n_ref> n_vec(n_ref);
70 scalar_seq_view<T_N_ref> N_vec(N_ref);
71 scalar_seq_view<T_theta_ref> theta_vec(theta_ref);
72 size_t size_theta = stan::math::size(theta);
73 size_t max_size_seq_view = max_size(n, N, theta);
74
76 for (size_t i = 0; i < size_theta; ++i) {
77 log1m_theta[i] = log1m(theta_vec.val(i));
78 }
79
81 for (size_t i = 0; i < max_size_seq_view; ++i) {
82 logp += binomial_coefficient_log(N_vec[i], n_vec[i]);
83 }
84 }
85
86 for (size_t i = 0; i < max_size_seq_view; ++i) {
87 if (N_vec[i] != 0) {
88 if (n_vec[i] == 0) {
89 logp += N_vec[i] * log1m_theta[i];
90 } else if (n_vec[i] == N_vec[i]) {
91 logp += n_vec[i] * log(theta_vec.val(i));
92 } else {
93 logp += n_vec[i] * log(theta_vec.val(i))
94 + (N_vec[i] - n_vec[i]) * log1m_theta[i];
95 }
96 }
97 }
98
100 if (size_theta == 1) {
101 T_partials_return sum_n = 0;
102 T_partials_return sum_N = 0;
103 for (size_t i = 0; i < max_size_seq_view; ++i) {
104 sum_n += n_vec[i];
105 sum_N += N_vec[i];
106 }
107 const T_partials_return theta_dbl = theta_vec.val(0);
108 if (sum_N != 0) {
109 if (sum_n == 0) {
110 partials<0>(ops_partials)[0] -= sum_N / (1.0 - theta_dbl);
111 } else if (sum_n == sum_N) {
112 partials<0>(ops_partials)[0] += sum_n / theta_dbl;
113 } else {
114 partials<0>(ops_partials)[0]
115 += sum_n / theta_dbl - (sum_N - sum_n) / (1.0 - theta_dbl);
116 }
117 }
118 } else {
119 for (size_t i = 0; i < max_size_seq_view; ++i) {
120 const T_partials_return theta_dbl = theta_vec.val(i);
121 if (N_vec[i] != 0) {
122 if (n_vec[i] == 0) {
123 partials<0>(ops_partials)[i] -= N_vec[i] / (1.0 - theta_dbl);
124 } else if (n_vec[i] == N_vec[i]) {
125 partials<0>(ops_partials)[i] += n_vec[i] / theta_dbl;
126 } else {
127 partials<0>(ops_partials)[i]
128 += n_vec[i] / theta_dbl
129 - (N_vec[i] - n_vec[i]) / (1.0 - theta_dbl);
130 }
131 }
132 }
133 }
134 }
135
136 return ops_partials.build(logp);
137}
138
139template <typename T_n, typename T_N, typename T_prob>
140inline return_type_t<T_prob> binomial_lpmf(const T_n& n, const T_N& N,
141 const T_prob& theta) {
142 return binomial_lpmf<false>(n, N, theta);
143}
144
145} // namespace math
146} // namespace stan
147#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_prob_cl > binomial_lpmf(const T_n_cl &n, const T_N_cl N, const T_prob_cl &theta)
Returns the log PMF for the binomial distribution evaluated at the specified 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
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
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
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.
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
fvar< T > log1m(const fvar< T > &x)
Definition log1m.hpp:12
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
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...