438 lines · plain
1// (C) Copyright Nick Thompson 2018.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_UNIVARIATE_STATISTICS_HPP7#define BOOST_MATH_TOOLS_UNIVARIATE_STATISTICS_HPP8 9#include <algorithm>10#include <iterator>11#include <tuple>12#include <boost/math/tools/assert.hpp>13#include <boost/math/tools/header_deprecated.hpp>14 15#include <boost/math/tools/is_standalone.hpp>16#ifndef BOOST_MATH_STANDALONE17#include <boost/config.hpp>18#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR19#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."20#endif21#endif22 23BOOST_MATH_HEADER_DEPRECATED("<boost/math/statistics/univariate_statistics.hpp>");24 25namespace boost::math::tools {26 27template<class ForwardIterator>28auto mean(ForwardIterator first, ForwardIterator last)29{30 using Real = typename std::iterator_traits<ForwardIterator>::value_type;31 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");32 if constexpr (std::is_integral<Real>::value)33 {34 double mu = 0;35 double i = 1;36 for(auto it = first; it != last; ++it) {37 mu = mu + (*it - mu)/i;38 i += 1;39 }40 return mu;41 }42 else if constexpr (std::is_same_v<typename std::iterator_traits<ForwardIterator>::iterator_category, std::random_access_iterator_tag>)43 {44 size_t elements = std::distance(first, last);45 Real mu0 = 0;46 Real mu1 = 0;47 Real mu2 = 0;48 Real mu3 = 0;49 Real i = 1;50 auto end = last - (elements % 4);51 for(auto it = first; it != end; it += 4) {52 Real inv = Real(1)/i;53 Real tmp0 = (*it - mu0);54 Real tmp1 = (*(it+1) - mu1);55 Real tmp2 = (*(it+2) - mu2);56 Real tmp3 = (*(it+3) - mu3);57 // please generate a vectorized fma here58 mu0 += tmp0*inv;59 mu1 += tmp1*inv;60 mu2 += tmp2*inv;61 mu3 += tmp3*inv;62 i += 1;63 }64 Real num1 = Real(elements - (elements %4))/Real(4);65 Real num2 = num1 + Real(elements % 4);66 67 for (auto it = end; it != last; ++it)68 {69 mu3 += (*it-mu3)/i;70 i += 1;71 }72 73 return (num1*(mu0+mu1+mu2) + num2*mu3)/Real(elements);74 }75 else76 {77 auto it = first;78 Real mu = *it;79 Real i = 2;80 while(++it != last)81 {82 mu += (*it - mu)/i;83 i += 1;84 }85 return mu;86 }87}88 89template<class Container>90inline auto mean(Container const & v)91{92 return mean(v.cbegin(), v.cend());93}94 95template<class ForwardIterator>96auto variance(ForwardIterator first, ForwardIterator last)97{98 using Real = typename std::iterator_traits<ForwardIterator>::value_type;99 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");100 // Higham, Accuracy and Stability, equation 1.6a and 1.6b:101 if constexpr (std::is_integral<Real>::value)102 {103 double M = *first;104 double Q = 0;105 double k = 2;106 for (auto it = std::next(first); it != last; ++it)107 {108 double tmp = *it - M;109 Q = Q + ((k-1)*tmp*tmp)/k;110 M = M + tmp/k;111 k += 1;112 }113 return Q/(k-1);114 }115 else116 {117 Real M = *first;118 Real Q = 0;119 Real k = 2;120 for (auto it = std::next(first); it != last; ++it)121 {122 Real tmp = (*it - M)/k;123 Q += k*(k-1)*tmp*tmp;124 M += tmp;125 k += 1;126 }127 return Q/(k-1);128 }129}130 131template<class Container>132inline auto variance(Container const & v)133{134 return variance(v.cbegin(), v.cend());135}136 137template<class ForwardIterator>138auto sample_variance(ForwardIterator first, ForwardIterator last)139{140 size_t n = std::distance(first, last);141 BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");142 return n*variance(first, last)/(n-1);143}144 145template<class Container>146inline auto sample_variance(Container const & v)147{148 return sample_variance(v.cbegin(), v.cend());149}150 151 152// Follows equation 1.5 of:153// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf154template<class ForwardIterator>155auto skewness(ForwardIterator first, ForwardIterator last)156{157 using Real = typename std::iterator_traits<ForwardIterator>::value_type;158 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute skewness.");159 if constexpr (std::is_integral<Real>::value)160 {161 double M1 = *first;162 double M2 = 0;163 double M3 = 0;164 double n = 2;165 for (auto it = std::next(first); it != last; ++it)166 {167 double delta21 = *it - M1;168 double tmp = delta21/n;169 M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);170 M2 = M2 + tmp*(n-1)*delta21;171 M1 = M1 + tmp;172 n += 1;173 }174 175 double var = M2/(n-1);176 if (var == 0)177 {178 // The limit is technically undefined, but the interpretation here is clear:179 // A constant dataset has no skewness.180 return static_cast<double>(0);181 }182 double skew = M3/(M2*sqrt(var));183 return skew;184 }185 else186 {187 Real M1 = *first;188 Real M2 = 0;189 Real M3 = 0;190 Real n = 2;191 for (auto it = std::next(first); it != last; ++it)192 {193 Real delta21 = *it - M1;194 Real tmp = delta21/n;195 M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);196 M2 += tmp*(n-1)*delta21;197 M1 += tmp;198 n += 1;199 }200 201 Real var = M2/(n-1);202 if (var == 0)203 {204 // The limit is technically undefined, but the interpretation here is clear:205 // A constant dataset has no skewness.206 return Real(0);207 }208 Real skew = M3/(M2*sqrt(var));209 return skew;210 }211}212 213template<class Container>214inline auto skewness(Container const & v)215{216 return skewness(v.cbegin(), v.cend());217}218 219// Follows equation 1.5/1.6 of:220// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf221template<class ForwardIterator>222auto first_four_moments(ForwardIterator first, ForwardIterator last)223{224 using Real = typename std::iterator_traits<ForwardIterator>::value_type;225 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the first four moments.");226 if constexpr (std::is_integral<Real>::value)227 {228 double M1 = *first;229 double M2 = 0;230 double M3 = 0;231 double M4 = 0;232 double n = 2;233 for (auto it = std::next(first); it != last; ++it)234 {235 double delta21 = *it - M1;236 double tmp = delta21/n;237 M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);238 M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);239 M2 = M2 + tmp*(n-1)*delta21;240 M1 = M1 + tmp;241 n += 1;242 }243 244 return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));245 }246 else247 {248 Real M1 = *first;249 Real M2 = 0;250 Real M3 = 0;251 Real M4 = 0;252 Real n = 2;253 for (auto it = std::next(first); it != last; ++it)254 {255 Real delta21 = *it - M1;256 Real tmp = delta21/n;257 M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);258 M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);259 M2 = M2 + tmp*(n-1)*delta21;260 M1 = M1 + tmp;261 n += 1;262 }263 264 return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));265 }266}267 268template<class Container>269inline auto first_four_moments(Container const & v)270{271 return first_four_moments(v.cbegin(), v.cend());272}273 274 275// Follows equation 1.6 of:276// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf277template<class ForwardIterator>278auto kurtosis(ForwardIterator first, ForwardIterator last)279{280 auto [M1, M2, M3, M4] = first_four_moments(first, last);281 if (M2 == 0)282 {283 return M2;284 }285 return M4/(M2*M2);286}287 288template<class Container>289inline auto kurtosis(Container const & v)290{291 return kurtosis(v.cbegin(), v.cend());292}293 294template<class ForwardIterator>295auto excess_kurtosis(ForwardIterator first, ForwardIterator last)296{297 return kurtosis(first, last) - 3;298}299 300template<class Container>301inline auto excess_kurtosis(Container const & v)302{303 return excess_kurtosis(v.cbegin(), v.cend());304}305 306 307template<class RandomAccessIterator>308auto median(RandomAccessIterator first, RandomAccessIterator last)309{310 size_t num_elems = std::distance(first, last);311 BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");312 if (num_elems & 1)313 {314 auto middle = first + (num_elems - 1)/2;315 std::nth_element(first, middle, last);316 return *middle;317 }318 else319 {320 auto middle = first + num_elems/2 - 1;321 std::nth_element(first, middle, last);322 std::nth_element(middle, middle+1, last);323 return (*middle + *(middle+1))/2;324 }325}326 327 328template<class RandomAccessContainer>329inline auto median(RandomAccessContainer & v)330{331 return median(v.begin(), v.end());332}333 334template<class RandomAccessIterator>335auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)336{337 using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;338 BOOST_MATH_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");339 340 std::sort(first, last);341 if constexpr (std::is_integral<Real>::value)342 {343 double i = 1;344 double num = 0;345 double denom = 0;346 for (auto it = first; it != last; ++it)347 {348 num += *it*i;349 denom += *it;350 ++i;351 }352 353 // If the l1 norm is zero, all elements are zero, so every element is the same.354 if (denom == 0)355 {356 return static_cast<double>(0);357 }358 359 return ((2*num)/denom - i)/(i-1);360 }361 else362 {363 Real i = 1;364 Real num = 0;365 Real denom = 0;366 for (auto it = first; it != last; ++it)367 {368 num += *it*i;369 denom += *it;370 ++i;371 }372 373 // If the l1 norm is zero, all elements are zero, so every element is the same.374 if (denom == 0)375 {376 return Real(0);377 }378 379 return ((2*num)/denom - i)/(i-1);380 }381}382 383template<class RandomAccessContainer>384inline auto gini_coefficient(RandomAccessContainer & v)385{386 return gini_coefficient(v.begin(), v.end());387}388 389template<class RandomAccessIterator>390inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)391{392 size_t n = std::distance(first, last);393 return n*gini_coefficient(first, last)/(n-1);394}395 396template<class RandomAccessContainer>397inline auto sample_gini_coefficient(RandomAccessContainer & v)398{399 return sample_gini_coefficient(v.begin(), v.end());400}401 402template<class RandomAccessIterator>403auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())404{405 using std::abs;406 using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;407 using std::isnan;408 if (isnan(center))409 {410 center = boost::math::tools::median(first, last);411 }412 size_t num_elems = std::distance(first, last);413 BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");414 auto comparator = [¢er](Real a, Real b) { return abs(a-center) < abs(b-center);};415 if (num_elems & 1)416 {417 auto middle = first + (num_elems - 1)/2;418 std::nth_element(first, middle, last, comparator);419 return abs(*middle);420 }421 else422 {423 auto middle = first + num_elems/2 - 1;424 std::nth_element(first, middle, last, comparator);425 std::nth_element(middle, middle+1, last, comparator);426 return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<Real>(2));427 }428}429 430template<class RandomAccessContainer>431inline auto median_absolute_deviation(RandomAccessContainer & v, typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())432{433 return median_absolute_deviation(v.begin(), v.end(), center);434}435 436}437#endif438