From c7acd2f391b850e34225051b63881fa0f9c4852d Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 15 Aug 2026 13:23:01 +0800 Subject: [PATCH 1/8] Use nth_element for quantile --- stan/math/prim/fun/quantile.hpp | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index a427cb6e9e9..e7effca6078 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -41,21 +41,25 @@ inline double quantile(const T& samples_vec, const double p) { Eigen::VectorXd x = as_array_or_scalar(samples_vec); check_not_nan("quantile", "samples_vec", x); - if (n_sample == 1) + if (n_sample == 1) { return x.coeff(0); - else if (p == 0.) - return *std::min_element(x.data(), x.data() + n_sample); - else if (p == 1.) - return *std::max_element(x.data(), x.data() + n_sample); - - double index = (n_sample - 1) * p; - size_t lo = std::floor(index); - size_t hi = std::ceil(index); + } else if (p == 0.) { + return x.minCoeff(); + } else if (p == 1.) { + return x.maxCoeff(); + } + + const double index = (n_sample - 1) * p; + const size_t lo = std::floor(index); - std::sort(x.data(), x.data() + n_sample, std::less()); + std::nth_element(x.data(), x.data() + lo, x.data() + n_sample); - double h = index - lo; - return (1 - h) * x.coeff(lo) + h * x.coeff(hi); + const double h = index - lo; + if (h == 0) { + return x.coeff(lo); + } + Eigen::Map hi_map(x.data() + lo + 1, n_sample - (lo + 1)); + return (1 - h) * x.coeff(lo) + h * hi_map.minCoeff(); } /** From a6beabf6ec117d240871c784976e981954044e68 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 15 Aug 2026 01:43:22 -0400 Subject: [PATCH 2/8] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/quantile.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index e7effca6078..d8ef18a0e28 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -48,7 +48,7 @@ inline double quantile(const T& samples_vec, const double p) { } else if (p == 1.) { return x.maxCoeff(); } - + const double index = (n_sample - 1) * p; const size_t lo = std::floor(index); @@ -58,7 +58,8 @@ inline double quantile(const T& samples_vec, const double p) { if (h == 0) { return x.coeff(lo); } - Eigen::Map hi_map(x.data() + lo + 1, n_sample - (lo + 1)); + Eigen::Map hi_map(x.data() + lo + 1, + n_sample - (lo + 1)); return (1 - h) * x.coeff(lo) + h * hi_map.minCoeff(); } From f893ed95e396a61a98fabd22f2cb867012841abc Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 15 Aug 2026 15:37:04 +0800 Subject: [PATCH 3/8] Update vectorised impl to take Eigen vecs --- stan/math/prim/fun/quantile.hpp | 27 +++++++++-------------- test/unit/math/prim/fun/quantile_test.cpp | 12 ++++++++-- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index d8ef18a0e28..37bdf232d08 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -1,6 +1,7 @@ #ifndef STAN_MATH_PRIM_FUN_QUANTILE_HPP #define STAN_MATH_PRIM_FUN_QUANTILE_HPP +#include #include #include #include @@ -82,8 +83,8 @@ inline double quantile(const T& samples_vec, const double p) { */ template * = nullptr, require_vector_vt* = nullptr, - require_std_vector_vt* = nullptr> -inline std::vector quantile(const T& samples_vec, const Tp& ps) { + require_vector_vt* = nullptr> +inline plain_type_t quantile(const T& samples_vec, const Tp& ps) { check_not_nan("quantile", "ps", ps); check_bounded("quantile", "ps", ps, 0, 1); @@ -96,25 +97,17 @@ inline std::vector quantile(const T& samples_vec, const Tp& ps) { Eigen::VectorXd x = as_array_or_scalar(samples_vec); check_not_nan("quantile", "samples_vec", x); - const auto& p = as_array_or_scalar(ps); - std::vector ret(n_ps, 0.0); + plain_type_t ret(n_ps); std::sort(x.data(), x.data() + n_sample, std::less()); - Eigen::ArrayXd index = (n_sample - 1) * p; + const Eigen::ArrayXd index = (n_sample - 1) * as_array_or_scalar(ps); + const Eigen::ArrayXd lo = index.floor(); + const Eigen::ArrayXd h = index - lo; for (size_t i = 0; i < n_ps; ++i) { - if (p[i] == 0.) { - ret[i] = x.coeff(0); - } else if (p[i] == 1.) { - ret[i] = x.coeff(n_sample - 1); - } else { - size_t lo = std::floor(index[i]); - size_t hi = std::ceil(index[i]); - - double h = index[i] - lo; - - ret[i] = (1 - h) * x.coeff(lo) + h * x.coeff(hi); - } + const size_t l = static_cast(lo.coeff(i)); + const size_t u = l + (h.coeff(i) > 0); + ret[i] = (1 - h.coeff(i)) * x.coeff(l) + h.coeff(i) * x.coeff(u); } return ret; } diff --git a/test/unit/math/prim/fun/quantile_test.cpp b/test/unit/math/prim/fun/quantile_test.cpp index a939a1098ac..6d94d286554 100644 --- a/test/unit/math/prim/fun/quantile_test.cpp +++ b/test/unit/math/prim/fun/quantile_test.cpp @@ -105,7 +105,7 @@ inline void test_quantile_double() { p[2] = 0.2; p[3] = 1; - std::vector ret = quantile(v, p); + T ret = quantile(v, p); EXPECT_FLOAT_EQ(ret[0], -0.28); EXPECT_FLOAT_EQ(ret[1], -0.196); EXPECT_FLOAT_EQ(ret[2], -0.112); @@ -145,7 +145,7 @@ inline void test_quantile_double() { // check size 1 first argument works T v1(1); v1[0] = -0.07; - std::vector ret1 = quantile(v1, p); + T ret1 = quantile(v1, p); EXPECT_FLOAT_EQ(ret1[0], -0.07); EXPECT_FLOAT_EQ(ret1[1], -0.07); EXPECT_FLOAT_EQ(ret1[2], -0.07); @@ -160,6 +160,14 @@ TEST(MathFunctions, quantileEigenVectorXdStdVecDouble) { test_quantile_double>(); } +TEST(MathFunctions, quantileEigenVectorXdEigenVectorXd) { + test_quantile_double(); +} + TEST(MathFunctions, quantileEigenRowVectorXdStdVecDouble) { test_quantile_double>(); } + +TEST(MathFunctions, quantileEigenRowVectorXdEigenRowVectorXd) { + test_quantile_double(); +} \ No newline at end of file From 12ff13e05036ae649f5fb64ae844e5178c6e57a6 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 15 Aug 2026 15:40:55 +0800 Subject: [PATCH 4/8] Remove stray includes --- stan/math/prim/fun/quantile.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index 37bdf232d08..d1c9fcba902 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -1,13 +1,11 @@ #ifndef STAN_MATH_PRIM_FUN_QUANTILE_HPP #define STAN_MATH_PRIM_FUN_QUANTILE_HPP -#include #include #include #include #include #include -#include namespace stan { namespace math { From 0b0123f3ddafb2766fc9e97103f60411d0784208 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 15 Aug 2026 15:48:30 +0800 Subject: [PATCH 5/8] cpplint --- test/unit/math/prim/fun/quantile_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/quantile_test.cpp b/test/unit/math/prim/fun/quantile_test.cpp index 6d94d286554..4b6308eb53b 100644 --- a/test/unit/math/prim/fun/quantile_test.cpp +++ b/test/unit/math/prim/fun/quantile_test.cpp @@ -170,4 +170,4 @@ TEST(MathFunctions, quantileEigenRowVectorXdStdVecDouble) { TEST(MathFunctions, quantileEigenRowVectorXdEigenRowVectorXd) { test_quantile_double(); -} \ No newline at end of file +} From 5c546e301d80077a10fb098e544eff8aa90ea285 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 18 Aug 2026 22:21:48 +0800 Subject: [PATCH 6/8] Must return same type as probs vec --- stan/math/prim/fun/quantile.hpp | 41 ++++++++++++----------- test/unit/math/prim/fun/quantile_test.cpp | 4 +-- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index d1c9fcba902..90d63250188 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -48,8 +48,9 @@ inline double quantile(const T& samples_vec, const double p) { return x.maxCoeff(); } - const double index = (n_sample - 1) * p; - const size_t lo = std::floor(index); + const size_t nm1 = (n_sample - 1); + const double index = nm1 * p; + const size_t lo = static_cast(index); std::nth_element(x.data(), x.data() + lo, x.data() + n_sample); @@ -57,9 +58,7 @@ inline double quantile(const T& samples_vec, const double p) { if (h == 0) { return x.coeff(lo); } - Eigen::Map hi_map(x.data() + lo + 1, - n_sample - (lo + 1)); - return (1 - h) * x.coeff(lo) + h * hi_map.minCoeff(); + return (1 - h) * x.coeff(lo) + h * x.tail(nm1 - lo).minCoeff(); } /** @@ -79,10 +78,12 @@ inline double quantile(const T& samples_vec, const double p) { * @throw std::invalid_argument If any of the values are NaN or size 0. * @throw std::domain_error If `p<0` or `p>1` for any p in ps. */ -template * = nullptr, +template , + require_all_vector_t* = nullptr, require_vector_vt* = nullptr, require_vector_vt* = nullptr> -inline plain_type_t quantile(const T& samples_vec, const Tp& ps) { +inline ReturnT quantile(const T& samples_vec, const Tp& ps) { check_not_nan("quantile", "ps", ps); check_bounded("quantile", "ps", ps, 0, 1); @@ -92,21 +93,23 @@ inline plain_type_t quantile(const T& samples_vec, const Tp& ps) { return {}; } - Eigen::VectorXd x = as_array_or_scalar(samples_vec); - check_not_nan("quantile", "samples_vec", x); + check_not_nan("quantile", "samples_vec", samples_vec); + Eigen::ArrayXd x = as_array_or_scalar(samples_vec); + std::sort(x.begin(), x.end()); - plain_type_t ret(n_ps); + ReturnT ret(n_ps); + Eigen::Map ret_map(ret.data(), n_ps); + Eigen::Map ps_map(ps.data(), n_ps); - std::sort(x.data(), x.data() + n_sample, std::less()); - const Eigen::ArrayXd index = (n_sample - 1) * as_array_or_scalar(ps); - const Eigen::ArrayXd lo = index.floor(); - const Eigen::ArrayXd h = index - lo; + const size_t nm1 = (n_sample - 1); + ret_map = ps_map.unaryExpr([&x, nm1](const double ps_i) { + const double idx_xi = nm1 * ps_i; + const size_t lo = static_cast(idx_xi); + const double h = idx_xi - lo; + + return (1 - h) * x[lo] + h * x[lo + (h != 0.0)]; + }); - for (size_t i = 0; i < n_ps; ++i) { - const size_t l = static_cast(lo.coeff(i)); - const size_t u = l + (h.coeff(i) > 0); - ret[i] = (1 - h.coeff(i)) * x.coeff(l) + h.coeff(i) * x.coeff(u); - } return ret; } diff --git a/test/unit/math/prim/fun/quantile_test.cpp b/test/unit/math/prim/fun/quantile_test.cpp index 4b6308eb53b..a2cb4dbc74c 100644 --- a/test/unit/math/prim/fun/quantile_test.cpp +++ b/test/unit/math/prim/fun/quantile_test.cpp @@ -105,7 +105,7 @@ inline void test_quantile_double() { p[2] = 0.2; p[3] = 1; - T ret = quantile(v, p); + Tp ret = quantile(v, p); EXPECT_FLOAT_EQ(ret[0], -0.28); EXPECT_FLOAT_EQ(ret[1], -0.196); EXPECT_FLOAT_EQ(ret[2], -0.112); @@ -145,7 +145,7 @@ inline void test_quantile_double() { // check size 1 first argument works T v1(1); v1[0] = -0.07; - T ret1 = quantile(v1, p); + Tp ret1 = quantile(v1, p); EXPECT_FLOAT_EQ(ret1[0], -0.07); EXPECT_FLOAT_EQ(ret1[1], -0.07); EXPECT_FLOAT_EQ(ret1[2], -0.07); From 9ff28ac21a10ecef7537474a83b673d16f9db5e9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 18 Aug 2026 23:34:06 +0800 Subject: [PATCH 7/8] Simplify vectorised version further --- stan/math/prim/fun/quantile.hpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index 90d63250188..81a4787ea54 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -94,21 +94,22 @@ inline ReturnT quantile(const T& samples_vec, const Tp& ps) { } check_not_nan("quantile", "samples_vec", samples_vec); - Eigen::ArrayXd x = as_array_or_scalar(samples_vec); - std::sort(x.begin(), x.end()); + plain_type_t x = samples_vec; + std::sort(x.begin(), x.end()); ReturnT ret(n_ps); - Eigen::Map ret_map(ret.data(), n_ps); - Eigen::Map ps_map(ps.data(), n_ps); const size_t nm1 = (n_sample - 1); - ret_map = ps_map.unaryExpr([&x, nm1](const double ps_i) { - const double idx_xi = nm1 * ps_i; - const size_t lo = static_cast(idx_xi); - const double h = idx_xi - lo; - return (1 - h) * x[lo] + h * x[lo + (h != 0.0)]; - }); + for (size_t i = 0; i < n_ps; i++) { + const double sample_xi = nm1 * ps[i]; + const size_t smpl_xi_int = static_cast(sample_xi); + const double smpl_xi_frc = sample_xi - smpl_xi_int; + ret[i] = (1 - smpl_xi_frc) * x[smpl_xi_int]; + if (smpl_xi_frc != 0.0) { + ret[i] += smpl_xi_frc * x[smpl_xi_int + 1]; + } + } return ret; } From e8cd75bb9da3486ba6389ab8bfefe141a2bc14f5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 19 Aug 2026 14:56:58 +0800 Subject: [PATCH 8/8] Fix expression evals --- stan/math/prim/fun/quantile.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/quantile.hpp b/stan/math/prim/fun/quantile.hpp index 81a4787ea54..f343e495419 100644 --- a/stan/math/prim/fun/quantile.hpp +++ b/stan/math/prim/fun/quantile.hpp @@ -93,9 +93,9 @@ inline ReturnT quantile(const T& samples_vec, const Tp& ps) { return {}; } - check_not_nan("quantile", "samples_vec", samples_vec); - plain_type_t x = samples_vec; + check_not_nan("quantile", "samples_vec", x); + std::sort(x.begin(), x.end()); ReturnT ret(n_ps);