Skip to content

Commit 3cef4b6

Browse files
lingiumlingium
committed
add beta_neg_binomial_lccdf
1 parent 9c7c3ff commit 3cef4b6

3 files changed

Lines changed: 248 additions & 0 deletions

File tree

‎stan/math/prim/prob.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include <stan/math/prim/prob/beta_lccdf.hpp>
2626
#include <stan/math/prim/prob/beta_lcdf.hpp>
2727
#include <stan/math/prim/prob/beta_lpdf.hpp>
28+
#include <stan/math/prim/prob/beta_neg_binomial_lccdf.hpp>
2829
#include <stan/math/prim/prob/beta_neg_binomial_lpmf.hpp>
2930
#include <stan/math/prim/prob/beta_proportion_ccdf_log.hpp>
3031
#include <stan/math/prim/prob/beta_proportion_cdf_log.hpp>

‎stan/math/prim/prob/beta_neg_binomial_lccdf.hpp‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
#ifndef STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_LCCDF_HPP
2+
#define STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_LCCDF_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/constants.hpp>
7+
#include <stan/math/prim/fun/digamma.hpp>
8+
#include <stan/math/prim/fun/hypergeometric_3F2.hpp>
9+
#include <stan/math/prim/fun/grad_F32.hpp>
10+
#include <stan/math/prim/fun/lbeta.hpp>
11+
#include <stan/math/prim/fun/lgamma.hpp>
12+
#include <stan/math/prim/fun/max_size.hpp>
13+
#include <stan/math/prim/fun/scalar_seq_view.hpp>
14+
#include <stan/math/prim/fun/size.hpp>
15+
#include <stan/math/prim/fun/size_zero.hpp>
16+
#include <stan/math/prim/functor/partials_propagator.hpp>
17+
#include <cmath>
18+
19+
namespace stan {
20+
namespace math {
21+
22+
/** \ingroup prob_dists
23+
* Returns the log CCDF of the Beta-Negative Binomial distribution with given
24+
* number of successes, prior success, and prior failure parameters.
25+
* Given containers of matching sizes, returns the log sum of probabilities.
26+
*
27+
* @tparam T_n type of failure parameter
28+
* @tparam T_r type of number of successes parameter
29+
* @tparam T_alpha type of prior success parameter
30+
* @tparam T_beta type of prior failure parameter
31+
*
32+
* @param n failure parameter
33+
* @param r Number of successes parameter
34+
* @param alpha prior success parameter
35+
* @param beta prior failure parameter
36+
* @return log probability or log sum of probabilities
37+
* @throw std::domain_error if r, alpha, or beta fails to be positive
38+
* @throw std::invalid_argument if container sizes mismatch
39+
*/
40+
template <typename T_n, typename T_r, typename T_alpha, typename T_beta>
41+
inline return_type_t<T_r, T_alpha, T_beta> beta_neg_binomial_lccdf(
42+
const T_n& n, const T_r& r, const T_alpha& alpha, const T_beta& beta) {
43+
using std::exp;
44+
using std::log;
45+
using T_partials_return = partials_return_t<T_n, T_r, T_alpha, T_beta>;
46+
using T_r_ref = ref_type_t<T_r>;
47+
using T_alpha_ref = ref_type_t<T_alpha>;
48+
using T_beta_ref = ref_type_t<T_beta>;
49+
static constexpr const char* function = "beta_neg_binomial_lccdf";
50+
check_consistent_sizes(
51+
function, "Failures variable", n, "Number of successes parameter", r,
52+
"Prior success parameter", alpha, "Prior failure parameter", beta);
53+
if (size_zero(n, r, alpha, beta)) {
54+
return 0;
55+
}
56+
57+
T_r_ref r_ref = r;
58+
T_alpha_ref alpha_ref = alpha;
59+
T_beta_ref beta_ref = beta;
60+
check_positive_finite(function, "Number of successes parameter", r_ref);
61+
check_positive_finite(function, "Prior success parameter", alpha_ref);
62+
check_positive_finite(function, "Prior failure parameter", beta_ref);
63+
64+
T_partials_return P(0.0);
65+
auto ops_partials = make_partials_propagator(r_ref, alpha_ref, beta_ref);
66+
67+
scalar_seq_view<T_n> n_vec(n);
68+
scalar_seq_view<T_r_ref> r_vec(r_ref);
69+
scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
70+
scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
71+
size_t size_n = stan::math::size(n);
72+
size_t max_size_seq_view = max_size(n, r, alpha, beta);
73+
74+
// Explicit return for extreme values
75+
// The gradients are technically ill-defined, but treated as zero
76+
for (size_t i = 0; i < size_n; i++) {
77+
if (n_vec.val(i) < 0) {
78+
return ops_partials.build(0.0);
79+
}
80+
}
81+
82+
for (size_t i = 0; i < max_size_seq_view; i++) {
83+
// Explicit return for extreme values
84+
// The gradients are technically ill-defined, but treated as zero
85+
if (n_vec.val(i) == std::numeric_limits<int>::max()) {
86+
return ops_partials.build(negative_infinity());
87+
}
88+
T_partials_return n_dbl = n_vec.val(i);
89+
T_partials_return r_dbl = r_vec.val(i);
90+
T_partials_return alpha_dbl = alpha_vec.val(i);
91+
T_partials_return beta_dbl = beta_vec.val(i);
92+
T_partials_return b_plus_n = beta_dbl + n_dbl;
93+
T_partials_return r_plus_n = r_dbl + n_dbl;
94+
T_partials_return a_plus_r = alpha_dbl + r_dbl;
95+
T_partials_return one = 1;
96+
T_partials_return precision = 1e-8; // default -6, set -8 to pass all tests
97+
98+
T_partials_return F
99+
= hypergeometric_3F2({one, b_plus_n + 1, r_plus_n + 1},
100+
{n_dbl + 2, a_plus_r + b_plus_n + 1}, one);
101+
T_partials_return C = lgamma(r_plus_n + 1) + lbeta(a_plus_r, b_plus_n + 1)
102+
- lgamma(r_dbl) - lbeta(alpha_dbl, beta_dbl)
103+
- lgamma(n_dbl + 2);
104+
T_partials_return ccdf = exp(C) * F;
105+
T_partials_return P_i = log(ccdf);
106+
P += P_i;
107+
108+
if (!is_constant_all<T_r, T_alpha, T_beta>::value) {
109+
T_partials_return digamma_n_r_alpha_beta
110+
= digamma(a_plus_r + b_plus_n + 1);
111+
T_partials_return dF[6];
112+
grad_F32(dF, one, b_plus_n + 1, r_plus_n + 1, n_dbl + 2,
113+
a_plus_r + b_plus_n + 1, one, precision, 1e5);
114+
115+
if constexpr (!is_constant<T_r>::value || !is_constant<T_alpha>::value) {
116+
T_partials_return digamma_r_alpha = digamma(a_plus_r);
117+
if constexpr (!is_constant_all<T_r>::value) {
118+
partials<0>(ops_partials)[i]
119+
+= digamma(r_plus_n + 1)
120+
+ (digamma_r_alpha - digamma_n_r_alpha_beta)
121+
+ (dF[2] + dF[4]) / F - digamma(r_dbl);
122+
}
123+
if constexpr (!is_constant_all<T_alpha>::value) {
124+
partials<1>(ops_partials)[i] += digamma_r_alpha
125+
- digamma_n_r_alpha_beta + dF[4] / F
126+
- digamma(alpha_dbl);
127+
}
128+
}
129+
130+
if constexpr (!is_constant<T_alpha>::value
131+
|| !is_constant<T_beta>::value) {
132+
T_partials_return digamma_alpha_beta = digamma(alpha_dbl + beta_dbl);
133+
if constexpr (!is_constant<T_alpha>::value) {
134+
partials<1>(ops_partials)[i] += digamma_alpha_beta;
135+
}
136+
if constexpr (!is_constant<T_beta>::value) {
137+
partials<2>(ops_partials)[i]
138+
+= digamma(b_plus_n + 1) - digamma_n_r_alpha_beta
139+
+ (dF[1] + dF[4]) / F
140+
- (digamma(beta_dbl) - digamma_alpha_beta);
141+
}
142+
}
143+
}
144+
}
145+
146+
return ops_partials.build(P);
147+
}
148+
149+
} // namespace math
150+
} // namespace stan
151+
#endif

‎test/prob/beta_neg_binomial/beta_neg_binomial_ccdf_log_test.hpp‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
// Arguments: Ints, Doubles, Doubles, Doubles
2+
#include <stan/math/prim/prob/beta_neg_binomial_lccdf.hpp>
3+
#include <stan/math/prim/fun/lbeta.hpp>
4+
#include <stan/math/prim/fun/lgamma.hpp>
5+
6+
using stan::math::var;
7+
using std::numeric_limits;
8+
using std::vector;
9+
10+
class AgradCcdfLogBetaNegBinomial : public AgradCcdfLogTest {
11+
public:
12+
void valid_values(vector<vector<double>>& parameters,
13+
vector<double>& ccdf_log) {
14+
vector<double> param(4);
15+
16+
param[0] = 10; // n
17+
param[1] = 5.5; // r
18+
param[2] = 2.5; // alpha
19+
param[3] = 0.5; // beta
20+
parameters.push_back(param);
21+
ccdf_log.push_back(std::log(1.0 - 0.967906252841089)); // expected ccdf_log
22+
}
23+
24+
void invalid_values(vector<size_t>& index, vector<double>& value) {
25+
// n
26+
27+
// r
28+
index.push_back(1U);
29+
value.push_back(0.0);
30+
31+
index.push_back(1U);
32+
value.push_back(-1.0);
33+
34+
index.push_back(1U);
35+
value.push_back(std::numeric_limits<double>::infinity());
36+
37+
// alpha
38+
index.push_back(2U);
39+
value.push_back(0.0);
40+
41+
index.push_back(2U);
42+
value.push_back(-1.0);
43+
44+
index.push_back(2U);
45+
value.push_back(std::numeric_limits<double>::infinity());
46+
47+
// beta
48+
index.push_back(3U);
49+
value.push_back(0.0);
50+
51+
index.push_back(3U);
52+
value.push_back(-1.0);
53+
54+
index.push_back(3U);
55+
value.push_back(std::numeric_limits<double>::infinity());
56+
}
57+
58+
// BOUND INCLUDED IN ORDER FOR TEST TO PASS WITH CURRENT FRAMEWORK
59+
bool has_lower_bound() { return false; }
60+
61+
bool has_upper_bound() { return false; }
62+
63+
template <typename T_n, typename T_r, typename T_size1, typename T_size2,
64+
typename T4, typename T5>
65+
stan::return_type_t<T_r, T_size1, T_size2> ccdf_log(const T_n& n,
66+
const T_r& r,
67+
const T_size1& alpha,
68+
const T_size2& beta,
69+
const T4&, const T5&) {
70+
return stan::math::beta_neg_binomial_lccdf(n, r, alpha, beta);
71+
}
72+
73+
template <typename T_n, typename T_r, typename T_size1, typename T_size2,
74+
typename T4, typename T5>
75+
stan::return_type_t<T_r, T_size1, T_size2> ccdf_log_function(
76+
const T_n& n, const T_r& r, const T_size1& alpha, const T_size2& beta,
77+
const T4&, const T5&) {
78+
using stan::math::lbeta;
79+
using stan::math::lgamma;
80+
using stan::math::log1m;
81+
using stan::math::log_sum_exp;
82+
using std::vector;
83+
84+
vector<stan::return_type_t<T_r, T_size1, T_size2>> lpmf_values;
85+
86+
for (int i = 0; i <= n; i++) {
87+
auto lpmf = lbeta(i + r, alpha + beta) - lbeta(r, alpha)
88+
+ lgamma(i + beta) - lgamma(i + 1) - lgamma(beta);
89+
lpmf_values.push_back(lpmf);
90+
}
91+
92+
auto log_cdf = log_sum_exp(lpmf_values);
93+
94+
return log1m(exp(log_cdf));
95+
}
96+
};

0 commit comments

Comments
 (0)