1185 lines · plain
1// (C) Copyright Nick Thompson 2018.2// (C) Copyright Matt Borland 2020.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP8#define BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP9 10#include <boost/math/statistics/detail/single_pass.hpp>11#include <boost/math/tools/config.hpp>12#include <boost/math/tools/assert.hpp>13#include <algorithm>14#include <iterator>15#include <tuple>16#include <cmath>17#include <vector>18#include <type_traits>19#include <utility>20#include <numeric>21#include <list>22 23#ifdef BOOST_MATH_EXEC_COMPATIBLE24#include <execution>25 26namespace boost::math::statistics {27 28template<class ExecutionPolicy, class ForwardIterator>29inline auto mean(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)30{31 using Real = typename std::iterator_traits<ForwardIterator>::value_type;32 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");33 34 if constexpr (std::is_integral_v<Real>)35 {36 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)37 {38 return detail::mean_sequential_impl<double>(first, last);39 }40 else41 {42 return std::reduce(exec, first, last, 0.0) / std::distance(first, last);43 }44 }45 else46 {47 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)48 {49 return detail::mean_sequential_impl<Real>(first, last);50 }51 else52 {53 return std::reduce(exec, first, last, Real(0.0)) / Real(std::distance(first, last));54 }55 }56}57 58template<class ExecutionPolicy, class Container>59inline auto mean(ExecutionPolicy&& exec, Container const & v)60{61 return mean(exec, std::cbegin(v), std::cend(v));62}63 64template<class ForwardIterator>65inline auto mean(ForwardIterator first, ForwardIterator last)66{67 return mean(std::execution::seq, first, last);68}69 70template<class Container>71inline auto mean(Container const & v)72{73 return mean(std::execution::seq, std::cbegin(v), std::cend(v));74}75 76template<class ExecutionPolicy, class ForwardIterator>77inline auto variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)78{79 using Real = typename std::iterator_traits<ForwardIterator>::value_type;80 81 if constexpr (std::is_integral_v<Real>)82 {83 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)84 {85 return std::get<2>(detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last));86 }87 else88 {89 const auto results = detail::first_four_moments_parallel_impl<std::tuple<double, double, double, double, double>>(first, last);90 return std::get<1>(results) / std::get<4>(results);91 }92 }93 else94 {95 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)96 {97 return std::get<2>(detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last));98 }99 else100 {101 const auto results = detail::first_four_moments_parallel_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last);102 return std::get<1>(results) / std::get<4>(results);103 }104 }105}106 107template<class ExecutionPolicy, class Container>108inline auto variance(ExecutionPolicy&& exec, Container const & v)109{110 return variance(exec, std::cbegin(v), std::cend(v));111}112 113template<class ForwardIterator>114inline auto variance(ForwardIterator first, ForwardIterator last)115{116 return variance(std::execution::seq, first, last);117}118 119template<class Container>120inline auto variance(Container const & v)121{122 return variance(std::execution::seq, std::cbegin(v), std::cend(v));123}124 125template<class ExecutionPolicy, class ForwardIterator>126inline auto sample_variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)127{128 const auto n = std::distance(first, last);129 BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");130 return n*variance(exec, first, last)/(n-1);131}132 133template<class ExecutionPolicy, class Container>134inline auto sample_variance(ExecutionPolicy&& exec, Container const & v)135{136 return sample_variance(exec, std::cbegin(v), std::cend(v));137}138 139template<class ForwardIterator>140inline auto sample_variance(ForwardIterator first, ForwardIterator last)141{142 return sample_variance(std::execution::seq, first, last);143}144 145template<class Container>146inline auto sample_variance(Container const & v)147{148 return sample_variance(std::execution::seq, std::cbegin(v), std::cend(v));149}150 151template<class ExecutionPolicy, class ForwardIterator>152inline auto mean_and_sample_variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)153{154 using Real = typename std::iterator_traits<ForwardIterator>::value_type;155 156 if constexpr (std::is_integral_v<Real>)157 {158 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)159 {160 const auto results = detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last);161 return std::make_pair(std::get<0>(results), std::get<2>(results)*std::get<3>(results)/(std::get<3>(results)-1.0));162 }163 else164 {165 const auto results = detail::first_four_moments_parallel_impl<std::tuple<double, double, double, double, double>>(first, last);166 return std::make_pair(std::get<0>(results), std::get<1>(results) / (std::get<4>(results)-1.0));167 }168 }169 else170 {171 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)172 {173 const auto results = detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last);174 return std::make_pair(std::get<0>(results), std::get<2>(results)*std::get<3>(results)/(std::get<3>(results)-Real(1)));175 }176 else177 {178 const auto results = detail::first_four_moments_parallel_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last);179 return std::make_pair(std::get<0>(results), std::get<1>(results) / (std::get<4>(results)-Real(1)));180 }181 }182}183 184template<class ExecutionPolicy, class Container>185inline auto mean_and_sample_variance(ExecutionPolicy&& exec, Container const & v)186{187 return mean_and_sample_variance(exec, std::cbegin(v), std::cend(v));188}189 190template<class ForwardIterator>191inline auto mean_and_sample_variance(ForwardIterator first, ForwardIterator last)192{193 return mean_and_sample_variance(std::execution::seq, first, last);194}195 196template<class Container>197inline auto mean_and_sample_variance(Container const & v)198{199 return mean_and_sample_variance(std::execution::seq, std::cbegin(v), std::cend(v));200}201 202template<class ExecutionPolicy, class ForwardIterator>203inline auto first_four_moments(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)204{205 using Real = typename std::iterator_traits<ForwardIterator>::value_type;206 207 if constexpr (std::is_integral_v<Real>)208 {209 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)210 {211 const auto results = detail::first_four_moments_sequential_impl<std::tuple<double, double, double, double, double>>(first, last);212 return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results),213 std::get<3>(results) / std::get<4>(results));214 }215 else216 {217 const auto results = detail::first_four_moments_parallel_impl<std::tuple<double, double, double, double, double>>(first, last);218 return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results),219 std::get<3>(results) / std::get<4>(results));220 }221 }222 else223 {224 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)225 {226 const auto results = detail::first_four_moments_sequential_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last);227 return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results),228 std::get<3>(results) / std::get<4>(results));229 }230 else231 {232 const auto results = detail::first_four_moments_parallel_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last);233 return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results),234 std::get<3>(results) / std::get<4>(results));235 }236 }237}238 239template<class ExecutionPolicy, class Container>240inline auto first_four_moments(ExecutionPolicy&& exec, Container const & v)241{242 return first_four_moments(exec, std::cbegin(v), std::cend(v));243}244 245template<class ForwardIterator>246inline auto first_four_moments(ForwardIterator first, ForwardIterator last)247{248 return first_four_moments(std::execution::seq, first, last);249}250 251template<class Container>252inline auto first_four_moments(Container const & v)253{254 return first_four_moments(std::execution::seq, std::cbegin(v), std::cend(v));255}256 257// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf258template<class ExecutionPolicy, class ForwardIterator>259inline auto skewness(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)260{261 using Real = typename std::iterator_traits<ForwardIterator>::value_type;262 using std::sqrt;263 264 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)265 {266 if constexpr (std::is_integral_v<Real>)267 {268 return detail::skewness_sequential_impl<double>(first, last);269 }270 else271 {272 return detail::skewness_sequential_impl<Real>(first, last);273 }274 }275 else276 {277 const auto [M1, M2, M3, M4] = first_four_moments(exec, first, last);278 const auto n = std::distance(first, last);279 const auto var = M2/(n-1);280 281 if (M2 == 0)282 {283 // The limit is technically undefined, but the interpretation here is clear:284 // A constant dataset has no skewness.285 if constexpr (std::is_integral_v<Real>)286 {287 return static_cast<double>(0);288 }289 else290 {291 return Real(0);292 }293 }294 else295 {296 return M3/(M2*sqrt(var)) / Real(2);297 }298 }299}300 301template<class ExecutionPolicy, class Container>302inline auto skewness(ExecutionPolicy&& exec, Container & v)303{304 return skewness(exec, std::cbegin(v), std::cend(v));305}306 307template<class ForwardIterator>308inline auto skewness(ForwardIterator first, ForwardIterator last)309{310 return skewness(std::execution::seq, first, last);311}312 313template<class Container>314inline auto skewness(Container const & v)315{316 return skewness(std::execution::seq, std::cbegin(v), std::cend(v));317}318 319// Follows equation 1.6 of:320// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf321template<class ExecutionPolicy, class ForwardIterator>322inline auto kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)323{324 const auto [M1, M2, M3, M4] = first_four_moments(exec, first, last);325 if (M2 == 0)326 {327 return M2;328 }329 return M4/(M2*M2);330}331 332template<class ExecutionPolicy, class Container>333inline auto kurtosis(ExecutionPolicy&& exec, Container const & v)334{335 return kurtosis(exec, std::cbegin(v), std::cend(v));336}337 338template<class ForwardIterator>339inline auto kurtosis(ForwardIterator first, ForwardIterator last)340{341 return kurtosis(std::execution::seq, first, last);342}343 344template<class Container>345inline auto kurtosis(Container const & v)346{347 return kurtosis(std::execution::seq, std::cbegin(v), std::cend(v));348}349 350template<class ExecutionPolicy, class ForwardIterator>351inline auto excess_kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)352{353 return kurtosis(exec, first, last) - 3;354}355 356template<class ExecutionPolicy, class Container>357inline auto excess_kurtosis(ExecutionPolicy&& exec, Container const & v)358{359 return excess_kurtosis(exec, std::cbegin(v), std::cend(v));360}361 362template<class ForwardIterator>363inline auto excess_kurtosis(ForwardIterator first, ForwardIterator last)364{365 return excess_kurtosis(std::execution::seq, first, last);366}367 368template<class Container>369inline auto excess_kurtosis(Container const & v)370{371 return excess_kurtosis(std::execution::seq, std::cbegin(v), std::cend(v));372}373 374 375template<class ExecutionPolicy, class RandomAccessIterator>376auto median(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)377{378 const auto num_elems = std::distance(first, last);379 BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");380 if (num_elems & 1)381 {382 auto middle = first + (num_elems - 1)/2;383 std::nth_element(exec, first, middle, last);384 return *middle;385 }386 else387 {388 auto middle = first + num_elems/2 - 1;389 std::nth_element(exec, first, middle, last);390 std::nth_element(exec, middle, middle+1, last);391 return (*middle + *(middle+1))/2;392 }393}394 395 396template<class ExecutionPolicy, class RandomAccessContainer>397inline auto median(ExecutionPolicy&& exec, RandomAccessContainer & v)398{399 return median(exec, std::begin(v), std::end(v));400}401 402template<class RandomAccessIterator>403inline auto median(RandomAccessIterator first, RandomAccessIterator last)404{405 return median(std::execution::seq, first, last);406}407 408template<class RandomAccessContainer>409inline auto median(RandomAccessContainer & v)410{411 return median(std::execution::seq, std::begin(v), std::end(v));412}413 414#if 0415//416// Parallel gini calculation is curently broken, see:417// https://github.com/boostorg/math/issues/585418// We will fix this at a later date, for now just use a serial implementation:419//420template<class ExecutionPolicy, class RandomAccessIterator>421inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)422{423 using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;424 425 if(!std::is_sorted(exec, first, last))426 {427 std::sort(exec, first, last);428 }429 430 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)431 {432 if constexpr (std::is_integral_v<Real>)433 {434 return detail::gini_coefficient_sequential_impl<double>(first, last);435 }436 else437 {438 return detail::gini_coefficient_sequential_impl<Real>(first, last);439 }440 }441 442 else if constexpr (std::is_integral_v<Real>)443 {444 return detail::gini_coefficient_parallel_impl<double>(exec, first, last);445 }446 447 else448 {449 return detail::gini_coefficient_parallel_impl<Real>(exec, first, last);450 }451}452#else453template<class ExecutionPolicy, class RandomAccessIterator>454inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)455{456 using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;457 458 if (!std::is_sorted(exec, first, last))459 {460 std::sort(exec, first, last);461 }462 463 if constexpr (std::is_integral_v<Real>)464 {465 return detail::gini_coefficient_sequential_impl<double>(first, last);466 }467 else468 {469 return detail::gini_coefficient_sequential_impl<Real>(first, last);470 }471}472#endif473 474template<class ExecutionPolicy, class RandomAccessContainer>475inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessContainer & v)476{477 return gini_coefficient(exec, std::begin(v), std::end(v));478}479 480template<class RandomAccessIterator>481inline auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)482{483 return gini_coefficient(std::execution::seq, first, last);484}485 486template<class RandomAccessContainer>487inline auto gini_coefficient(RandomAccessContainer & v)488{489 return gini_coefficient(std::execution::seq, std::begin(v), std::end(v));490}491 492template<class ExecutionPolicy, class RandomAccessIterator>493inline auto sample_gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)494{495 const auto n = std::distance(first, last);496 return n*gini_coefficient(exec, first, last)/(n-1);497}498 499template<class ExecutionPolicy, class RandomAccessContainer>500inline auto sample_gini_coefficient(ExecutionPolicy&& exec, RandomAccessContainer & v)501{502 return sample_gini_coefficient(exec, std::begin(v), std::end(v));503}504 505template<class RandomAccessIterator>506inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)507{508 return sample_gini_coefficient(std::execution::seq, first, last);509}510 511template<class RandomAccessContainer>512inline auto sample_gini_coefficient(RandomAccessContainer & v)513{514 return sample_gini_coefficient(std::execution::seq, std::begin(v), std::end(v));515}516 517template<class ExecutionPolicy, class RandomAccessIterator>518auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last,519 typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())520{521 using std::abs;522 using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;523 using std::isnan;524 if (isnan(center))525 {526 center = boost::math::statistics::median(exec, first, last);527 }528 const auto num_elems = std::distance(first, last);529 BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");530 auto comparator = [¢er](Real a, Real b) { return abs(a-center) < abs(b-center);};531 if (num_elems & 1)532 {533 auto middle = first + (num_elems - 1)/2;534 std::nth_element(exec, first, middle, last, comparator);535 return abs(*middle-center);536 }537 else538 {539 auto middle = first + num_elems/2 - 1;540 std::nth_element(exec, first, middle, last, comparator);541 std::nth_element(exec, middle, middle+1, last, comparator);542 return (abs(*middle-center) + abs(*(middle+1)-center))/abs(static_cast<Real>(2));543 }544}545 546template<class ExecutionPolicy, class RandomAccessContainer>547inline auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessContainer & v,548 typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())549{550 return median_absolute_deviation(exec, std::begin(v), std::end(v), center);551}552 553template<class RandomAccessIterator>554inline auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last,555 typename RandomAccessIterator::value_type center=std::numeric_limits<typename RandomAccessIterator::value_type>::quiet_NaN())556{557 return median_absolute_deviation(std::execution::seq, first, last, center);558}559 560template<class RandomAccessContainer>561inline auto median_absolute_deviation(RandomAccessContainer & v,562 typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())563{564 return median_absolute_deviation(std::execution::seq, std::begin(v), std::end(v), center);565}566 567template<class ExecutionPolicy, class ForwardIterator>568auto interquartile_range(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)569{570 using Real = typename std::iterator_traits<ForwardIterator>::value_type;571 static_assert(!std::is_integral_v<Real>, "Integer values have not yet been implemented.");572 auto m = std::distance(first,last);573 BOOST_MATH_ASSERT_MSG(m >= 3, "At least 3 samples are required to compute the interquartile range.");574 auto k = m/4;575 auto j = m - (4*k);576 // m = 4k+j.577 // If j = 0 or j = 1, then there are an even number of samples below the median, and an even number above the median.578 // Then we must average adjacent elements to get the quartiles.579 // If j = 2 or j = 3, there are an odd number of samples above and below the median, these elements may be directly extracted to get the quartiles.580 581 if (j==2 || j==3)582 {583 auto q1 = first + k;584 auto q3 = first + 3*k + j - 1;585 std::nth_element(exec, first, q1, last);586 Real Q1 = *q1;587 std::nth_element(exec, q1, q3, last);588 Real Q3 = *q3;589 return Q3 - Q1;590 } else {591 // j == 0 or j==1:592 auto q1 = first + k - 1;593 auto q3 = first + 3*k - 1 + j;594 std::nth_element(exec, first, q1, last);595 Real a = *q1;596 std::nth_element(exec, q1, q1 + 1, last);597 Real b = *(q1 + 1);598 Real Q1 = (a+b)/2;599 std::nth_element(exec, q1, q3, last);600 a = *q3;601 std::nth_element(exec, q3, q3 + 1, last);602 b = *(q3 + 1);603 Real Q3 = (a+b)/2;604 return Q3 - Q1;605 }606}607 608template<class ExecutionPolicy, class RandomAccessContainer>609inline auto interquartile_range(ExecutionPolicy&& exec, RandomAccessContainer & v)610{611 return interquartile_range(exec, std::begin(v), std::end(v));612}613 614template<class RandomAccessIterator>615inline auto interquartile_range(RandomAccessIterator first, RandomAccessIterator last)616{617 return interquartile_range(std::execution::seq, first, last);618}619 620template<class RandomAccessContainer>621inline auto interquartile_range(RandomAccessContainer & v)622{623 return interquartile_range(std::execution::seq, std::begin(v), std::end(v));624}625 626template<class ExecutionPolicy, class ForwardIterator, class OutputIterator>627inline OutputIterator mode(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last, OutputIterator output)628{629 if(!std::is_sorted(exec, first, last))630 {631 if constexpr (std::is_same_v<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>)632 {633 std::sort(exec, first, last);634 }635 else636 {637 BOOST_MATH_ASSERT("Data must be sorted for sequential mode calculation");638 }639 }640 641 return detail::mode_impl(first, last, output);642}643 644template<class ExecutionPolicy, class Container, class OutputIterator>645inline OutputIterator mode(ExecutionPolicy&& exec, Container & v, OutputIterator output)646{647 return mode(exec, std::begin(v), std::end(v), output);648}649 650template<class ForwardIterator, class OutputIterator>651inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output)652{653 return mode(std::execution::seq, first, last, output);654}655 656// Requires enable_if_t to not clash with impl that returns std::list657// Very ugly. std::is_execution_policy_v returns false for the std::execution objects and decltype of the objects (e.g. std::execution::seq)658template<class Container, class OutputIterator, std::enable_if_t<!std::is_convertible_v<std::execution::sequenced_policy, Container> &&659 !std::is_convertible_v<std::execution::parallel_unsequenced_policy, Container> &&660 !std::is_convertible_v<std::execution::parallel_policy, Container>661 #if __cpp_lib_execution > 201900662 && !std::is_convertible_v<std::execution::unsequenced_policy, Container>663 #endif664 , bool> = true>665inline OutputIterator mode(Container & v, OutputIterator output)666{667 return mode(std::execution::seq, std::begin(v), std::end(v), output);668}669 670// std::list is the return type for the proposed STL stats library671 672template<class ExecutionPolicy, class ForwardIterator, class Real = typename std::iterator_traits<ForwardIterator>::value_type>673inline auto mode(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)674{675 std::list<Real> modes;676 mode(exec, first, last, std::inserter(modes, modes.begin()));677 return modes;678}679 680template<class ExecutionPolicy, class Container>681inline auto mode(ExecutionPolicy&& exec, Container & v)682{683 return mode(exec, std::begin(v), std::end(v));684}685 686template<class ForwardIterator>687inline auto mode(ForwardIterator first, ForwardIterator last)688{689 return mode(std::execution::seq, first, last);690}691 692template<class Container>693inline auto mode(Container & v)694{695 return mode(std::execution::seq, std::begin(v), std::end(v));696}697 698} // Namespace boost::math::statistics699 700#else // Backwards compatible bindings for C++11 or execution is not implemented701 702namespace boost { namespace math { namespace statistics {703 704template<bool B, class T = void>705using enable_if_t = typename std::enable_if<B, T>::type;706 707template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,708 enable_if_t<std::is_integral<Real>::value, bool> = true>709inline double mean(const ForwardIterator first, const ForwardIterator last)710{711 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");712 return detail::mean_sequential_impl<double>(first, last);713}714 715template<class Container, typename Real = typename Container::value_type,716 enable_if_t<std::is_integral<Real>::value, bool> = true>717inline double mean(const Container& c)718{719 return mean(std::begin(c), std::end(c));720}721 722template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,723 enable_if_t<!std::is_integral<Real>::value, bool> = true>724inline Real mean(const ForwardIterator first, const ForwardIterator last)725{726 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");727 return detail::mean_sequential_impl<Real>(first, last);728}729 730template<class Container, typename Real = typename Container::value_type,731 enable_if_t<!std::is_integral<Real>::value, bool> = true>732inline Real mean(const Container& c)733{734 return mean(std::begin(c), std::end(c));735}736 737template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,738 enable_if_t<std::is_integral<Real>::value, bool> = true>739inline double variance(const ForwardIterator first, const ForwardIterator last)740{741 return std::get<2>(detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last));742}743 744template<class Container, typename Real = typename Container::value_type,745 enable_if_t<std::is_integral<Real>::value, bool> = true>746inline double variance(const Container& c)747{748 return variance(std::begin(c), std::end(c));749}750 751template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,752 enable_if_t<!std::is_integral<Real>::value, bool> = true>753inline Real variance(const ForwardIterator first, const ForwardIterator last)754{755 return std::get<2>(detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last));756 757}758 759template<class Container, typename Real = typename Container::value_type,760 enable_if_t<!std::is_integral<Real>::value, bool> = true>761inline Real variance(const Container& c)762{763 return variance(std::begin(c), std::end(c));764}765 766template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,767 enable_if_t<std::is_integral<Real>::value, bool> = true>768inline double sample_variance(const ForwardIterator first, const ForwardIterator last)769{770 const auto n = std::distance(first, last);771 BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");772 return n*variance(first, last)/(n-1);773}774 775template<class Container, typename Real = typename Container::value_type,776 enable_if_t<std::is_integral<Real>::value, bool> = true>777inline double sample_variance(const Container& c)778{779 return sample_variance(std::begin(c), std::end(c));780}781 782template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,783 enable_if_t<!std::is_integral<Real>::value, bool> = true>784inline Real sample_variance(const ForwardIterator first, const ForwardIterator last)785{786 const auto n = std::distance(first, last);787 BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");788 return n*variance(first, last)/(n-1);789}790 791template<class Container, typename Real = typename Container::value_type,792 enable_if_t<!std::is_integral<Real>::value, bool> = true>793inline Real sample_variance(const Container& c)794{795 return sample_variance(std::begin(c), std::end(c));796}797 798template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,799 enable_if_t<std::is_integral<Real>::value, bool> = true>800inline std::pair<double, double> mean_and_sample_variance(const ForwardIterator first, const ForwardIterator last)801{802 const auto results = detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last);803 return std::make_pair(std::get<0>(results), std::get<3>(results)*std::get<2>(results)/(std::get<3>(results)-1.0));804}805 806template<class Container, typename Real = typename Container::value_type,807 enable_if_t<std::is_integral<Real>::value, bool> = true>808inline std::pair<double, double> mean_and_sample_variance(const Container& c)809{810 return mean_and_sample_variance(std::begin(c), std::end(c));811}812 813template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,814 enable_if_t<!std::is_integral<Real>::value, bool> = true>815inline std::pair<Real, Real> mean_and_sample_variance(const ForwardIterator first, const ForwardIterator last)816{817 const auto results = detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last);818 return std::make_pair(std::get<0>(results), std::get<3>(results)*std::get<2>(results)/(std::get<3>(results)-Real(1)));819}820 821template<class Container, typename Real = typename Container::value_type,822 enable_if_t<!std::is_integral<Real>::value, bool> = true>823inline std::pair<Real, Real> mean_and_sample_variance(const Container& c)824{825 return mean_and_sample_variance(std::begin(c), std::end(c));826}827 828template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,829 enable_if_t<std::is_integral<Real>::value, bool> = true>830inline std::tuple<double, double, double, double> first_four_moments(const ForwardIterator first, const ForwardIterator last)831{832 const auto results = detail::first_four_moments_sequential_impl<std::tuple<double, double, double, double, double>>(first, last);833 return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results),834 std::get<3>(results) / std::get<4>(results));835}836 837template<class Container, typename Real = typename Container::value_type,838 enable_if_t<std::is_integral<Real>::value, bool> = true>839inline std::tuple<double, double, double, double> first_four_moments(const Container& c)840{841 return first_four_moments(std::begin(c), std::end(c));842}843 844template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,845 enable_if_t<!std::is_integral<Real>::value, bool> = true>846inline std::tuple<Real, Real, Real, Real> first_four_moments(const ForwardIterator first, const ForwardIterator last)847{848 const auto results = detail::first_four_moments_sequential_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last);849 return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results),850 std::get<3>(results) / std::get<4>(results));851}852 853template<class Container, typename Real = typename Container::value_type,854 enable_if_t<!std::is_integral<Real>::value, bool> = true>855inline std::tuple<Real, Real, Real, Real> first_four_moments(const Container& c)856{857 return first_four_moments(std::begin(c), std::end(c));858}859 860template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,861 enable_if_t<std::is_integral<Real>::value, bool> = true>862inline double skewness(const ForwardIterator first, const ForwardIterator last)863{864 return detail::skewness_sequential_impl<double>(first, last);865}866 867template<class Container, typename Real = typename Container::value_type,868 enable_if_t<std::is_integral<Real>::value, bool> = true>869inline double skewness(const Container& c)870{871 return skewness(std::begin(c), std::end(c));872}873 874template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,875 enable_if_t<!std::is_integral<Real>::value, bool> = true>876inline Real skewness(const ForwardIterator first, const ForwardIterator last)877{878 return detail::skewness_sequential_impl<Real>(first, last);879}880 881template<class Container, typename Real = typename Container::value_type,882 enable_if_t<!std::is_integral<Real>::value, bool> = true>883inline Real skewness(const Container& c)884{885 return skewness(std::begin(c), std::end(c));886}887 888template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,889 enable_if_t<std::is_integral<Real>::value, bool> = true>890inline double kurtosis(const ForwardIterator first, const ForwardIterator last)891{892 std::tuple<double, double, double, double> M = first_four_moments(first, last);893 894 if(std::get<1>(M) == 0)895 {896 return std::get<1>(M);897 }898 else899 {900 return std::get<3>(M)/(std::get<1>(M)*std::get<1>(M));901 }902}903 904template<class Container, typename Real = typename Container::value_type,905 enable_if_t<std::is_integral<Real>::value, bool> = true>906inline double kurtosis(const Container& c)907{908 return kurtosis(std::begin(c), std::end(c));909}910 911template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,912 enable_if_t<!std::is_integral<Real>::value, bool> = true>913inline Real kurtosis(const ForwardIterator first, const ForwardIterator last)914{915 std::tuple<Real, Real, Real, Real> M = first_four_moments(first, last);916 917 if(std::get<1>(M) == 0)918 {919 return std::get<1>(M);920 }921 else922 {923 return std::get<3>(M)/(std::get<1>(M)*std::get<1>(M));924 }925}926 927template<class Container, typename Real = typename Container::value_type,928 enable_if_t<!std::is_integral<Real>::value, bool> = true>929inline Real kurtosis(const Container& c)930{931 return kurtosis(std::begin(c), std::end(c));932}933 934template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,935 enable_if_t<std::is_integral<Real>::value, bool> = true>936inline double excess_kurtosis(const ForwardIterator first, const ForwardIterator last)937{938 return kurtosis(first, last) - 3;939}940 941template<class Container, typename Real = typename Container::value_type,942 enable_if_t<std::is_integral<Real>::value, bool> = true>943inline double excess_kurtosis(const Container& c)944{945 return excess_kurtosis(std::begin(c), std::end(c));946}947 948template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,949 enable_if_t<!std::is_integral<Real>::value, bool> = true>950inline Real excess_kurtosis(const ForwardIterator first, const ForwardIterator last)951{952 return kurtosis(first, last) - 3;953}954 955template<class Container, typename Real = typename Container::value_type,956 enable_if_t<!std::is_integral<Real>::value, bool> = true>957inline Real excess_kurtosis(const Container& c)958{959 return excess_kurtosis(std::begin(c), std::end(c));960}961 962template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type>963Real median(RandomAccessIterator first, RandomAccessIterator last)964{965 const auto num_elems = std::distance(first, last);966 BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");967 if (num_elems & 1)968 {969 auto middle = first + (num_elems - 1)/2;970 std::nth_element(first, middle, last);971 return *middle;972 }973 else974 {975 auto middle = first + num_elems/2 - 1;976 std::nth_element(first, middle, last);977 std::nth_element(middle, middle+1, last);978 return (*middle + *(middle+1))/2;979 }980}981 982template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type>983inline Real median(RandomAccessContainer& c)984{985 return median(std::begin(c), std::end(c));986}987 988template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,989 enable_if_t<std::is_integral<Real>::value, bool> = true>990inline double gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)991{992 if(!std::is_sorted(first, last))993 {994 std::sort(first, last);995 }996 997 return detail::gini_coefficient_sequential_impl<double>(first, last);998}999 1000template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,1001 enable_if_t<std::is_integral<Real>::value, bool> = true>1002inline double gini_coefficient(RandomAccessContainer& c)1003{1004 return gini_coefficient(std::begin(c), std::end(c));1005}1006 1007template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,1008 enable_if_t<!std::is_integral<Real>::value, bool> = true>1009inline Real gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)1010{1011 if(!std::is_sorted(first, last))1012 {1013 std::sort(first, last);1014 }1015 1016 return detail::gini_coefficient_sequential_impl<Real>(first, last);1017}1018 1019template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,1020 enable_if_t<!std::is_integral<Real>::value, bool> = true>1021inline Real gini_coefficient(RandomAccessContainer& c)1022{1023 return gini_coefficient(std::begin(c), std::end(c));1024}1025 1026template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,1027 enable_if_t<std::is_integral<Real>::value, bool> = true>1028inline double sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)1029{1030 const auto n = std::distance(first, last);1031 return n*gini_coefficient(first, last)/(n-1);1032}1033 1034template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,1035 enable_if_t<std::is_integral<Real>::value, bool> = true>1036inline double sample_gini_coefficient(RandomAccessContainer& c)1037{1038 return sample_gini_coefficient(std::begin(c), std::end(c));1039}1040 1041template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,1042 enable_if_t<!std::is_integral<Real>::value, bool> = true>1043inline Real sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)1044{1045 const auto n = std::distance(first, last);1046 return n*gini_coefficient(first, last)/(n-1);1047}1048 1049template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,1050 enable_if_t<!std::is_integral<Real>::value, bool> = true>1051inline Real sample_gini_coefficient(RandomAccessContainer& c)1052{1053 return sample_gini_coefficient(std::begin(c), std::end(c));1054}1055 1056template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type>1057Real median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last,1058 typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())1059{1060 using std::abs;1061 using std::isnan;1062 if (isnan(center))1063 {1064 center = boost::math::statistics::median(first, last);1065 }1066 const auto num_elems = std::distance(first, last);1067 BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");1068 auto comparator = [¢er](Real a, Real b) { return abs(a-center) < abs(b-center);};1069 if (num_elems & 1)1070 {1071 auto middle = first + (num_elems - 1)/2;1072 std::nth_element(first, middle, last, comparator);1073 return abs(*middle-center);1074 }1075 else1076 {1077 auto middle = first + num_elems/2 - 1;1078 std::nth_element(first, middle, last, comparator);1079 std::nth_element(middle, middle+1, last, comparator);1080 return (abs(*middle-center) + abs(*(middle+1)-center))/abs(static_cast<Real>(2));1081 }1082}1083 1084template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type>1085inline Real median_absolute_deviation(RandomAccessContainer& c,1086 typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())1087{1088 return median_absolute_deviation(std::begin(c), std::end(c), center);1089}1090 1091template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>1092Real interquartile_range(ForwardIterator first, ForwardIterator last)1093{1094 static_assert(!std::is_integral<Real>::value, "Integer values have not yet been implemented.");1095 auto m = std::distance(first,last);1096 BOOST_MATH_ASSERT_MSG(m >= 3, "At least 3 samples are required to compute the interquartile range.");1097 auto k = m/4;1098 auto j = m - (4*k);1099 // m = 4k+j.1100 // If j = 0 or j = 1, then there are an even number of samples below the median, and an even number above the median.1101 // Then we must average adjacent elements to get the quartiles.1102 // If j = 2 or j = 3, there are an odd number of samples above and below the median, these elements may be directly extracted to get the quartiles.1103 1104 if (j==2 || j==3)1105 {1106 auto q1 = first + k;1107 auto q3 = first + 3*k + j - 1;1108 std::nth_element(first, q1, last);1109 Real Q1 = *q1;1110 std::nth_element(q1, q3, last);1111 Real Q3 = *q3;1112 return Q3 - Q1;1113 }1114 else1115 {1116 // j == 0 or j==1:1117 auto q1 = first + k - 1;1118 auto q3 = first + 3*k - 1 + j;1119 std::nth_element(first, q1, last);1120 Real a = *q1;1121 std::nth_element(q1, q1 + 1, last);1122 Real b = *(q1 + 1);1123 Real Q1 = (a+b)/2;1124 std::nth_element(q1, q3, last);1125 a = *q3;1126 std::nth_element(q3, q3 + 1, last);1127 b = *(q3 + 1);1128 Real Q3 = (a+b)/2;1129 return Q3 - Q1;1130 }1131}1132 1133template<class Container, typename Real = typename Container::value_type>1134Real interquartile_range(Container& c)1135{1136 return interquartile_range(std::begin(c), std::end(c));1137}1138 1139template<class ForwardIterator, class OutputIterator,1140 enable_if_t<std::is_same<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>::value, bool> = true>1141inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output)1142{1143 if(!std::is_sorted(first, last))1144 {1145 std::sort(first, last);1146 }1147 1148 return detail::mode_impl(first, last, output);1149}1150 1151template<class ForwardIterator, class OutputIterator,1152 enable_if_t<!std::is_same<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>::value, bool> = true>1153inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output)1154{1155 if(!std::is_sorted(first, last))1156 {1157 BOOST_MATH_ASSERT("Data must be sorted for mode calculation");1158 }1159 1160 return detail::mode_impl(first, last, output);1161}1162 1163template<class Container, class OutputIterator>1164inline OutputIterator mode(Container& c, OutputIterator output)1165{1166 return mode(std::begin(c), std::end(c), output);1167}1168 1169template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>1170inline std::list<Real> mode(ForwardIterator first, ForwardIterator last)1171{1172 std::list<Real> modes;1173 mode(first, last, std::inserter(modes, modes.begin()));1174 return modes;1175}1176 1177template<class Container, typename Real = typename Container::value_type>1178inline std::list<Real> mode(Container& c)1179{1180 return mode(std::begin(c), std::end(c));1181}1182}}}1183#endif1184#endif // BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP1185