quantile.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_QUANTILE_HPP
#define STAN_MATH_PRIM_FUN_QUANTILE_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/prim/fun/as_array_or_scalar.hpp>
#include <algorithm>
 
namespace stan {
namespace math {
 
template <typename T, require_vector_t<T>* = nullptr,
          require_vector_vt<std::is_arithmetic, T>* = nullptr>
inline double quantile(const T& samples_vec, const double p) {
  check_not_nan("quantile", "p", p);
  check_bounded("quantile", "p", p, 0, 1);
 
  const size_t n_sample = samples_vec.size();
  if (n_sample == 0) {
    return {};
  }
 
  Eigen::VectorXd x = as_array_or_scalar(samples_vec);
  check_not_nan("quantile", "samples_vec", x);
 
  if (n_sample == 1) {
    return x.coeff(0);
  } else if (p == 0.) {
    return x.minCoeff();
  } else if (p == 1.) {
    return x.maxCoeff();
  }
 
  const size_t nm1 = (n_sample - 1);
  const double index = nm1 * p;
  const size_t lo = static_cast<size_t>(index);
 
  std::nth_element(x.data(), x.data() + lo, x.data() + n_sample);
 
  const double h = index - lo;
  if (h == 0) {
    return x.coeff(lo);
  }
  return (1 - h) * x.coeff(lo) + h * x.tail(nm1 - lo).minCoeff();
}
 
template <typename T, typename Tp,
          typename ReturnT = promote_scalar_t<double, Tp>,
          require_all_vector_t<T, Tp>* = nullptr,
          require_vector_vt<std::is_arithmetic, T>* = nullptr,
          require_vector_vt<std::is_arithmetic, Tp>* = nullptr>
inline ReturnT quantile(const T& samples_vec, const Tp& ps) {
  check_not_nan("quantile", "ps", ps);
  check_bounded("quantile", "ps", ps, 0, 1);
 
  const size_t n_sample = samples_vec.size();
  const size_t n_ps = ps.size();
  if (n_ps == 0 || n_sample == 0) {
    return {};
  }
 
  plain_type_t<T> x = samples_vec;
  check_not_nan("quantile", "samples_vec", x);
 
  std::sort(x.begin(), x.end());
  ReturnT ret(n_ps);
 
  const size_t nm1 = (n_sample - 1);
 
  for (size_t i = 0; i < n_ps; i++) {
    const double sample_xi = nm1 * ps[i];
    const size_t smpl_xi_int = static_cast<size_t>(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;
}
 
}  // namespace math
}  // namespace stan
 
#endif