math/distributions/empirical_cumulative_distribution_function.hpp

0001 //  Copyright Nick Thompson 2019.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_DISTRIBUTIONS_EMPIRICAL_CUMULATIVE_DISTRIBUTION_FUNCTION_HPP
0007 #define BOOST_MATH_DISTRIBUTIONS_EMPIRICAL_CUMULATIVE_DISTRIBUTION_FUNCTION_HPP
0008 #include <algorithm>
0009 #include <iterator>
0010 #include <stdexcept>
0011 #include <type_traits>
0012 #include <utility>
0013
0014 #include <boost/math/tools/is_standalone.hpp>
0015 #ifndef BOOST_MATH_STANDALONE
0016 #include <boost/config.hpp>
0017 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0018 #error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
0019 #endif
0020 #endif
0021
0022 namespace boost { namespace math{
0023
0024 template<class RandomAccessContainer>
0025 class empirical_cumulative_distribution_function {
0026     using Real = typename RandomAccessContainer::value_type;
0027 public:
0028     empirical_cumulative_distribution_function(RandomAccessContainer && v, bool sorted = false)
0029     {
0030         if (v.size() == 0) {
0031             throw std::domain_error("At least one sample is required to compute an empirical CDF.");
0032         }
0033         m_v = std::move(v);
0034         if (!sorted) {
0035             std::sort(m_v.begin(), m_v.end());
0036         }
0037     }
0038
0039     auto operator()(Real x) const {
0040        if constexpr (std::is_integral_v<Real>)
0041        {
0042          if (x < m_v[0]) {
0043            return static_cast<double>(0);
0044          }
0045          if (x >= m_v[m_v.size()-1]) {
0046            return static_cast<double>(1);
0047          }
0048          auto it = std::upper_bound(m_v.begin(), m_v.end(), x);
0049          return static_cast<double>(std::distance(m_v.begin(), it))/static_cast<double>(m_v.size());
0050        }
0051        else
0052        {
0053          if (x < m_v[0]) {
0054            return Real(0);
0055          }
0056          if (x >= m_v[m_v.size()-1]) {
0057            return Real(1);
0058          }
0059          auto it = std::upper_bound(m_v.begin(), m_v.end(), x);
0060          return static_cast<Real>(std::distance(m_v.begin(), it))/static_cast<Real>(m_v.size());
0061       }
0062     }
0063
0064     RandomAccessContainer&& return_data() {
0065         return std::move(m_v);
0066     }
0067
0068 private:
0069     RandomAccessContainer m_v;
0070 };
0071
0072 }}
0073 #endif