Automatic Differentiation
 
Loading...
Searching...
No Matches
quantile.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_QUANTILE_HPP
2#define STAN_MATH_PRIM_FUN_QUANTILE_HPP
3
8#include <algorithm>
9#include <vector>
10
11namespace stan {
12namespace math {
13
30template <typename T, require_vector_t<T>* = nullptr,
31 require_vector_vt<std::is_arithmetic, T>* = nullptr>
32inline double quantile(const T& samples_vec, const double p) {
33 check_not_nan("quantile", "p", p);
34 check_bounded("quantile", "p", p, 0, 1);
35
36 const size_t n_sample = samples_vec.size();
37 if (n_sample == 0) {
38 return {};
39 }
40
41 Eigen::VectorXd x = as_array_or_scalar(samples_vec);
42 check_not_nan("quantile", "samples_vec", x);
43
44 if (n_sample == 1)
45 return x.coeff(0);
46 else if (p == 0.)
47 return *std::min_element(x.data(), x.data() + n_sample);
48 else if (p == 1.)
49 return *std::max_element(x.data(), x.data() + n_sample);
50
51 double index = (n_sample - 1) * p;
52 size_t lo = std::floor(index);
53 size_t hi = std::ceil(index);
54
55 std::sort(x.data(), x.data() + n_sample, std::less<double>());
56
57 double h = index - lo;
58 return (1 - h) * x.coeff(lo) + h * x.coeff(hi);
59}
60
78template <typename T, typename Tp, require_all_vector_t<T, Tp>* = nullptr,
79 require_vector_vt<std::is_arithmetic, T>* = nullptr,
80 require_std_vector_vt<std::is_arithmetic, Tp>* = nullptr>
81inline std::vector<double> quantile(const T& samples_vec, const Tp& ps) {
82 check_not_nan("quantile", "ps", ps);
83 check_bounded("quantile", "ps", ps, 0, 1);
84
85 const size_t n_sample = samples_vec.size();
86 const size_t n_ps = ps.size();
87 if (n_ps == 0 || n_sample == 0) {
88 return {};
89 }
90
91 Eigen::VectorXd x = as_array_or_scalar(samples_vec);
92 check_not_nan("quantile", "samples_vec", x);
93
94 const auto& p = as_array_or_scalar(ps);
95 std::vector<double> ret(n_ps, 0.0);
96
97 std::sort(x.data(), x.data() + n_sample, std::less<double>());
98 Eigen::ArrayXd index = (n_sample - 1) * p;
99
100 for (size_t i = 0; i < n_ps; ++i) {
101 if (p[i] == 0.) {
102 ret[i] = x.coeff(0);
103 } else if (p[i] == 1.) {
104 ret[i] = x.coeff(n_sample - 1);
105 } else {
106 size_t lo = std::floor(index[i]);
107 size_t hi = std::ceil(index[i]);
108
109 double h = index[i] - lo;
110
111 ret[i] = (1 - h) * x.coeff(lo) + h * x.coeff(hi);
112 }
113 }
114 return ret;
115}
116
117} // namespace math
118} // namespace stan
119
120#endif
T as_array_or_scalar(T &&v)
Returns specified input value.
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Check if the value is between the low and high values, inclusively.
double quantile(const T &samples_vec, const double p)
Return sample quantiles corresponding to the given probabilities.
Definition quantile.hpp:32
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...