Automatic Differentiation
 
Loading...
Searching...
No Matches
neg_binomial_lpmf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_NEG_BINOMIAL_LPMF_HPP
2#define STAN_MATH_PRIM_PROB_NEG_BINOMIAL_LPMF_HPP
3
18#include <cmath>
19
20namespace stan {
21namespace math {
22
23namespace internal {
24// Exposing to let me use in tests
25// The current tests fail for 1e8 and pass for 1e9, so setting to 1e10
26constexpr double neg_binomial_alpha_cutoff = 1e10;
27} // namespace internal
28
29// NegBinomial(n|alpha, beta) [alpha > 0; beta > 0; n >= 0]
30template <bool propto, typename T_n, typename T_shape, typename T_inv_scale,
32 T_n, T_shape, T_inv_scale>* = nullptr>
34 const T_shape& alpha,
35 const T_inv_scale& beta) {
36 using T_partials_return = partials_return_t<T_n, T_shape, T_inv_scale>;
37 using std::log;
38 using T_n_ref = ref_type_t<T_n>;
39 using T_alpha_ref = ref_type_t<T_shape>;
40 using T_beta_ref = ref_type_t<T_inv_scale>;
41 static constexpr const char* function = "neg_binomial_lpmf";
42 check_consistent_sizes(function, "Failures variable", n, "Shape parameter",
43 alpha, "Inverse scale parameter", beta);
44 T_n_ref n_ref = n;
45 T_alpha_ref alpha_ref = alpha;
46 T_beta_ref beta_ref = beta;
47 check_nonnegative(function, "Failures variable", n_ref);
48 check_positive_finite(function, "Shape parameter", alpha_ref);
49 check_positive_finite(function, "Inverse scale parameter", beta_ref);
50
51 if (size_zero(n, alpha, beta)) {
52 return 0.0;
53 }
55 return 0.0;
56 }
57
58 T_partials_return logp(0.0);
59 auto ops_partials = make_partials_propagator(alpha_ref, beta_ref);
60
61 scalar_seq_view<T_n_ref> n_vec(n_ref);
62 scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
63 scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
64 size_t size_alpha = stan::math::size(alpha);
65 size_t size_beta = stan::math::size(beta);
66 size_t size_alpha_beta = max_size(alpha, beta);
67 size_t max_size_seq_view = max_size(n, alpha, beta);
68
69 VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
70 digamma_alpha(size_alpha);
72 for (size_t i = 0; i < size_alpha; ++i) {
73 digamma_alpha[i] = digamma(alpha_vec.val(i));
74 }
75 }
76
79 for (size_t i = 0; i < size_beta; ++i) {
80 const T_partials_return beta_dbl = beta_vec.val(i);
81 log1p_inv_beta[i] = log1p(inv(beta_dbl));
82 log1p_beta[i] = log1p(beta_dbl);
83 }
84
86 T_shape, T_inv_scale>
87 lambda_m_alpha_over_1p_beta(size_alpha_beta);
89 for (size_t i = 0; i < size_alpha_beta; ++i) {
90 const T_partials_return alpha_dbl = alpha_vec.val(i);
91 const T_partials_return beta_dbl = beta_vec.val(i);
92 lambda_m_alpha_over_1p_beta[i]
93 = alpha_dbl / beta_dbl - alpha_dbl / (1 + beta_dbl);
94 }
95 }
96
97 for (size_t i = 0; i < max_size_seq_view; i++) {
98 const T_partials_return alpha_dbl = alpha_vec.val(i);
99 const T_partials_return beta_dbl = beta_vec.val(i);
100
102 if (n_vec[i] != 0) {
103 logp += binomial_coefficient_log(n_vec[i] + alpha_dbl - 1.0,
104 alpha_dbl - 1.0);
105 }
106 }
107 logp -= alpha_dbl * log1p_inv_beta[i] + n_vec[i] * log1p_beta[i];
108
110 partials<0>(ops_partials)[i] += digamma(alpha_dbl + n_vec[i])
111 - digamma_alpha[i] - log1p_inv_beta[i];
112 }
114 partials<1>(ops_partials)[i]
115 += lambda_m_alpha_over_1p_beta[i] - n_vec[i] / (beta_dbl + 1.0);
116 }
117 }
118
119 return ops_partials.build(logp);
120}
121
122template <typename T_n, typename T_shape, typename T_inv_scale>
124 const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
125 return neg_binomial_lpmf<false>(n, alpha, beta);
126}
127
128} // namespace math
129} // namespace stan
130#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_shape_cl, T_inv_scale_cl > neg_binomial_lpmf(const T_n_cl &n, const T_shape_cl &alpha, const T_inv_scale_cl &beta)
The log of the negative binomial density for the specified scalars given the specified mean(s) and de...
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
constexpr double neg_binomial_alpha_cutoff
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_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
fvar< T > log1p(const fvar< T > &x)
Definition log1p.hpp:12
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
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...