From 949fddd0809db23209d1c13406e6ad8381ed29ff Mon Sep 17 00:00:00 2001 From: Zhi Ling <1336265834@qq.com> Date: Fri, 1 Nov 2024 16:39:04 +0800 Subject: [PATCH 1/2] add beta_neg_binomial_rng --- stan/math/prim/prob.hpp | 1 + stan/math/prim/prob/beta_neg_binomial_rng.hpp | 67 +++++++++++++++++++ .../math/prim/prob/beta_neg_binomial_test.cpp | 55 +++++++++++++++ 3 files changed, 123 insertions(+) create mode 100644 stan/math/prim/prob/beta_neg_binomial_rng.hpp create mode 100644 test/unit/math/prim/prob/beta_neg_binomial_test.cpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 919f63b62c0..e04c7e87172 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/beta_neg_binomial_rng.hpp b/stan/math/prim/prob/beta_neg_binomial_rng.hpp new file mode 100644 index 00000000000..7cd566f8640 --- /dev/null +++ b/stan/math/prim/prob/beta_neg_binomial_rng.hpp @@ -0,0 +1,67 @@ +#ifndef STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_RNG_HPP +#define STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_RNG_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return a beta-negative binomial random variate with given number of + * successes, prior success, and prior failure parameters, using the given + * random number generator. + * + * r, alpha, and beta can each be a scalar or a one-dimensional container. Any + * non-scalar inputs must be the same size. + * + * @tparam T_r type of number of successes parameter + * @tparam T_alpha type of prior success parameter + * @tparam T_beta type of prior failure parameter + * @tparam RNG type of random number generator + * + * @param r (Sequence of) number of successes parameter(s) + * @param alpha (Sequence of) positive success parameter(s) + * @param beta (Sequence of) positive failure parameter(s) + * @param rng random number generator + * @return (Sequence of) beta-binomial random variate(s) + * @throw std::domain_error if r is negative, or alpha or beta are nonpositive + * @throw std::invalid_argument if non-scalar arguments are of different sizes + */ +template +inline typename VectorBuilder::type +beta_neg_binomial_rng(const T_r &r, const T_alpha &alpha, const T_beta &beta, + RNG &rng) { + using T_r_ref = ref_type_t; + using T_alpha_ref = ref_type_t; + using T_beta_ref = ref_type_t; + static constexpr const char *function = "beta_neg_binomial_rng"; + check_consistent_sizes(function, "Number of successes parameter", r, + "Prior success parameter", alpha, + "Prior failure parameter", beta); + + T_r_ref r_ref = r; + T_alpha_ref alpha_ref = alpha; + T_beta_ref beta_ref = beta; + check_positive_finite(function, "Number of successes parameter", r_ref); + check_positive_finite(function, "Prior success parameter", alpha_ref); + check_positive_finite(function, "Prior failure parameter", beta_ref); + + using T_p = decltype(beta_rng(alpha_ref, beta_ref, rng)); + T_p p = beta_rng(alpha_ref, beta_ref, rng); + + scalar_seq_view p_vec(p); + size_t size_p = stan::math::size(p); + VectorBuilder odds_ratio_p(size_p); + for (size_t n = 0; n < size_p; ++n) { + odds_ratio_p[n] = stan::math::exp(stan::math::log(p_vec.val(n)) + - stan::math::log((1 - p_vec.val(n)))); + } + + return neg_binomial_rng(r_ref, odds_ratio_p.data(), rng); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/prim/prob/beta_neg_binomial_test.cpp b/test/unit/math/prim/prob/beta_neg_binomial_test.cpp new file mode 100644 index 00000000000..1a4a848daf3 --- /dev/null +++ b/test/unit/math/prim/prob/beta_neg_binomial_test.cpp @@ -0,0 +1,55 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +class BetaNegBinomialTestRig : public VectorIntRNGTestRig { + public: + BetaNegBinomialTestRig() + : VectorIntRNGTestRig(10000, 10, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, + {1.1, 3.1, 8.1}, {1, 3, 8}, {-3.0, -2.0, 0.0}, + {-3, -2, 0}, {3.1, 4.1, 5.1}, {3, 4, 5}, + {-2.1, -0.5, 0.0}, {-3, -1, 0}, {1.1, 3.1, 8.1}, + {1, 3, 8}, {-3.0, -2.0, 0.0}, {-3, -2, 0}) {} + + template + auto generate_samples(const T1& N, const T2& alpha, const T3& beta, + T_rng& rng) const { + return stan::math::beta_neg_binomial_rng(N, alpha, beta, rng); + } + + template + double pmf(int y, T1 r, double alpha, double beta) const { + return std::exp(stan::math::beta_neg_binomial_lpmf(y, r, alpha, beta)); + } +}; + +TEST(ProbDistributionsBetaNegBinomial, errorCheck) { + check_dist_throws_int_first_argument(BetaNegBinomialTestRig()); +} + +TEST(ProbDistributionsBetaNegBinomial, distributionCheck) { + check_counts_int_real_real(BetaNegBinomialTestRig()); +} + +TEST(ProbDistributionBetaNegBinomial, error_check) { + boost::random::mt19937 rng; + EXPECT_NO_THROW(stan::math::beta_neg_binomial_rng(4, 0.6, 2.0, rng)); + + EXPECT_THROW(stan::math::beta_neg_binomial_rng(-4, 0.6, 2, rng), + std::domain_error); + EXPECT_THROW(stan::math::beta_neg_binomial_rng(4, -0.6, 2, rng), + std::domain_error); + EXPECT_THROW(stan::math::beta_neg_binomial_rng(4, 0.6, -2, rng), + std::domain_error); + EXPECT_THROW(stan::math::beta_neg_binomial_rng( + 4, stan::math::positive_infinity(), 2, rng), + std::domain_error); + EXPECT_THROW(stan::math::beta_neg_binomial_rng( + 4, 0.6, stan::math::positive_infinity(), rng), + std::domain_error); +} From a225dd6731960f488f30d17141eafb5cc65f37f8 Mon Sep 17 00:00:00 2001 From: Zhi Ling <1336265834@qq.com> Date: Tue, 19 Nov 2024 15:21:19 +0800 Subject: [PATCH 2/2] fix issues --- stan/math/prim/prob/beta_neg_binomial_rng.hpp | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/prob/beta_neg_binomial_rng.hpp b/stan/math/prim/prob/beta_neg_binomial_rng.hpp index 7cd566f8640..19d0e49f7cf 100644 --- a/stan/math/prim/prob/beta_neg_binomial_rng.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_rng.hpp @@ -9,7 +9,7 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Return a beta-negative binomial random variate with given number of + * Return a beta-negative binomial random variate with the given number of * successes, prior success, and prior failure parameters, using the given * random number generator. * @@ -22,17 +22,16 @@ namespace math { * @tparam RNG type of random number generator * * @param r (Sequence of) number of successes parameter(s) - * @param alpha (Sequence of) positive success parameter(s) - * @param beta (Sequence of) positive failure parameter(s) + * @param alpha (Sequence of) prior success parameter(s) + * @param beta (Sequence of) prior failure parameter(s) * @param rng random number generator - * @return (Sequence of) beta-binomial random variate(s) - * @throw std::domain_error if r is negative, or alpha or beta are nonpositive + * @return (Sequence of) beta-negative binomial random variate(s) + * @throw std::domain_error if r, alpha, or beta are nonpositive * @throw std::invalid_argument if non-scalar arguments are of different sizes */ -template -inline typename VectorBuilder::type -beta_neg_binomial_rng(const T_r &r, const T_alpha &alpha, const T_beta &beta, - RNG &rng) { +template +inline auto beta_neg_binomial_rng(const T_r &r, const T_alpha &alpha, + const T_beta &beta, RNG &rng) { using T_r_ref = ref_type_t; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; @@ -55,8 +54,8 @@ beta_neg_binomial_rng(const T_r &r, const T_alpha &alpha, const T_beta &beta, size_t size_p = stan::math::size(p); VectorBuilder odds_ratio_p(size_p); for (size_t n = 0; n < size_p; ++n) { - odds_ratio_p[n] = stan::math::exp(stan::math::log(p_vec.val(n)) - - stan::math::log((1 - p_vec.val(n)))); + odds_ratio_p[n] + = stan::math::exp(stan::math::log(p_vec.val(n)) - log1m(p_vec.val(n))); } return neg_binomial_rng(r_ref, odds_ratio_p.data(), rng);