From a936073aa6f449026739726bc229bbedc49be6b4 Mon Sep 17 00:00:00 2001 From: Ling Zhi Date: Sun, 22 Feb 2026 18:56:36 +0800 Subject: [PATCH 1/4] add yule_simon_rng --- stan/math/prim/prob.hpp | 4 ++ stan/math/prim/prob/yule_simon_rng.hpp | 51 +++++++++++++++++ test/unit/math/prim/prob/yule_simon_test.cpp | 58 ++++++++++++++++++++ 3 files changed, 113 insertions(+) create mode 100755 stan/math/prim/prob/yule_simon_rng.hpp create mode 100644 test/unit/math/prim/prob/yule_simon_test.cpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 0bdff9f5f70..308b084a165 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -320,5 +320,9 @@ #include #include #include +#include +#include +#include #include +#include #endif diff --git a/stan/math/prim/prob/yule_simon_rng.hpp b/stan/math/prim/prob/yule_simon_rng.hpp new file mode 100755 index 00000000000..307660d4278 --- /dev/null +++ b/stan/math/prim/prob/yule_simon_rng.hpp @@ -0,0 +1,51 @@ +#ifndef STAN_MATH_PRIM_PROB_YULE_SIMON_RNG_HPP +#define STAN_MATH_PRIM_PROB_YULE_SIMON_RNG_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return a yule-simon random variate with the given shape parameter, + * using the given random number generator. + * + * alpha can be a scalar or a one-dimensional container. + * + * @tparam T_alpha type of shape parameter + * @tparam RNG type of random number generator + * + * @param alpha (Sequence of) shape parameter(s) + * @param rng random number generator + * @return (Sequence of) yule-simon random variate(s) + * @throw std::domain_error if alpha is nonpositive + */ +template +inline auto yule_simon_rng(const T_alpha &alpha, RNG &rng) { + using T_alpha_ref = ref_type_t; + static constexpr const char *function = "yule_simon_rng"; + + T_alpha_ref alpha_ref = alpha; + check_positive_finite(function, "Shape parameter", alpha_ref); + + using T_w = decltype(exponential_rng(alpha_ref, rng)); + T_w w = exponential_rng(alpha_ref, rng); + + scalar_seq_view w_vec(w); + size_t size_w = stan::math::size(w); + + VectorBuilder output(size_w); + for (size_t n = 0; n < size_w; ++n) { + double p = exp(-w_vec.val(n)); + double odds_ratio_p = exp(log(p) - log1m(p)); + output[n] = neg_binomial_rng(1.0, odds_ratio_p, rng) + 1; + } + + return output.data(); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/prim/prob/yule_simon_test.cpp b/test/unit/math/prim/prob/yule_simon_test.cpp new file mode 100644 index 00000000000..7223bf9845c --- /dev/null +++ b/test/unit/math/prim/prob/yule_simon_test.cpp @@ -0,0 +1,58 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class YuleSimonTestRig : public VectorIntRNGTestRig { + public: + YuleSimonTestRig() + : VectorIntRNGTestRig(10000, 10, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, + {2.1, 4.1, 10.1}, {1, 2, 3}, {-3.0, -2.0, 0.0}, + {-3, -1, 0}) {} + + template + auto generate_samples(const T1& alpha, const T2&, const T3&, + T_rng& rng) const { + return stan::math::yule_simon_rng(alpha, rng); + } + + template + double pmf(int y, T1 alpha, double, double) const { + return std::exp(stan::math::yule_simon_lpmf(y, alpha)); + } +}; + +TEST(ProbDistributionsYuleSimon, errorCheck) { + check_dist_throws_all_types(YuleSimonTestRig()); +} + +TEST(ProbDistributionsYuleSimon, distributionCheck) { + check_counts_real(YuleSimonTestRig()); +} + +TEST(ProbDistributionsYuleSimon, error_check) { + boost::random::mt19937 rng; + + // Valid parameters should not throw + EXPECT_NO_THROW(stan::math::yule_simon_rng(1.0, rng)); + EXPECT_NO_THROW(stan::math::yule_simon_rng(2.0, rng)); + EXPECT_NO_THROW(stan::math::yule_simon_rng(10.0, rng)); + + // Invalid parameters should throw domain_error + EXPECT_THROW(stan::math::yule_simon_rng(-0.5, rng), std::domain_error); + EXPECT_THROW(stan::math::yule_simon_rng(0.0, rng), std::domain_error); + EXPECT_THROW(stan::math::yule_simon_rng(-10.0, rng), std::domain_error); + + // Infinity should throw + EXPECT_THROW(stan::math::yule_simon_rng(stan::math::positive_infinity(), rng), + std::domain_error); + + // NaN should throw + EXPECT_THROW(stan::math::yule_simon_rng(stan::math::not_a_number(), rng), + std::domain_error); +} From 3bd34314ec30790786355e0bf6ed0b45aa30a0db Mon Sep 17 00:00:00 2001 From: Ling Zhi Date: Sun, 22 Feb 2026 20:25:18 +0800 Subject: [PATCH 2/4] fix test --- stan/math/prim/prob/yule_simon_rng.hpp | 6 ++++-- test/unit/math/prim/prob/yule_simon_test.cpp | 4 ---- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/prob/yule_simon_rng.hpp b/stan/math/prim/prob/yule_simon_rng.hpp index 307660d4278..5d8dcd72efb 100755 --- a/stan/math/prim/prob/yule_simon_rng.hpp +++ b/stan/math/prim/prob/yule_simon_rng.hpp @@ -2,6 +2,7 @@ #define STAN_MATH_PRIM_PROB_YULE_SIMON_RNG_HPP #include +#include #include #include @@ -38,8 +39,9 @@ inline auto yule_simon_rng(const T_alpha &alpha, RNG &rng) { VectorBuilder output(size_w); for (size_t n = 0; n < size_w; ++n) { - double p = exp(-w_vec.val(n)); - double odds_ratio_p = exp(log(p) - log1m(p)); + double p = stan::math::exp(-w_vec.val(n)); + double odds_ratio_p + = stan::math::exp(stan::math::log(p) - stan::math::log1m(p)); output[n] = neg_binomial_rng(1.0, odds_ratio_p, rng) + 1; } diff --git a/test/unit/math/prim/prob/yule_simon_test.cpp b/test/unit/math/prim/prob/yule_simon_test.cpp index 7223bf9845c..8935d088494 100644 --- a/test/unit/math/prim/prob/yule_simon_test.cpp +++ b/test/unit/math/prim/prob/yule_simon_test.cpp @@ -38,21 +38,17 @@ TEST(ProbDistributionsYuleSimon, distributionCheck) { TEST(ProbDistributionsYuleSimon, error_check) { boost::random::mt19937 rng; - // Valid parameters should not throw EXPECT_NO_THROW(stan::math::yule_simon_rng(1.0, rng)); EXPECT_NO_THROW(stan::math::yule_simon_rng(2.0, rng)); EXPECT_NO_THROW(stan::math::yule_simon_rng(10.0, rng)); - // Invalid parameters should throw domain_error EXPECT_THROW(stan::math::yule_simon_rng(-0.5, rng), std::domain_error); EXPECT_THROW(stan::math::yule_simon_rng(0.0, rng), std::domain_error); EXPECT_THROW(stan::math::yule_simon_rng(-10.0, rng), std::domain_error); - // Infinity should throw EXPECT_THROW(stan::math::yule_simon_rng(stan::math::positive_infinity(), rng), std::domain_error); - // NaN should throw EXPECT_THROW(stan::math::yule_simon_rng(stan::math::not_a_number(), rng), std::domain_error); } From 060a0df889dfab81f7599661561a71055e3d65f8 Mon Sep 17 00:00:00 2001 From: Ling Zhi Date: Sun, 22 Feb 2026 20:43:23 +0800 Subject: [PATCH 3/4] fix test again --- stan/math/prim/prob/yule_simon_rng.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/stan/math/prim/prob/yule_simon_rng.hpp b/stan/math/prim/prob/yule_simon_rng.hpp index 5d8dcd72efb..e51450c63c8 100755 --- a/stan/math/prim/prob/yule_simon_rng.hpp +++ b/stan/math/prim/prob/yule_simon_rng.hpp @@ -2,6 +2,8 @@ #define STAN_MATH_PRIM_PROB_YULE_SIMON_RNG_HPP #include +#include +#include #include #include #include From 405322bd0502ac523a86cafb1c3b10df540fb8cc Mon Sep 17 00:00:00 2001 From: Ling Zhi Date: Tue, 24 Feb 2026 18:01:52 +0800 Subject: [PATCH 4/4] fix issue --- stan/math/prim/prob/yule_simon_rng.hpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/prob/yule_simon_rng.hpp b/stan/math/prim/prob/yule_simon_rng.hpp index e51450c63c8..dd5e33bc7c6 100755 --- a/stan/math/prim/prob/yule_simon_rng.hpp +++ b/stan/math/prim/prob/yule_simon_rng.hpp @@ -33,16 +33,14 @@ inline auto yule_simon_rng(const T_alpha &alpha, RNG &rng) { T_alpha_ref alpha_ref = alpha; check_positive_finite(function, "Shape parameter", alpha_ref); - using T_w = decltype(exponential_rng(alpha_ref, rng)); - T_w w = exponential_rng(alpha_ref, rng); + auto w = exponential_rng(alpha_ref, rng); + scalar_seq_view w_vec(w); - scalar_seq_view w_vec(w); size_t size_w = stan::math::size(w); - VectorBuilder output(size_w); for (size_t n = 0; n < size_w; ++n) { - double p = stan::math::exp(-w_vec.val(n)); - double odds_ratio_p + const double p = stan::math::exp(-w_vec[n]); + const double odds_ratio_p = stan::math::exp(stan::math::log(p) - stan::math::log1m(p)); output[n] = neg_binomial_rng(1.0, odds_ratio_p, rng) + 1; }