402 lines · plain
1// (C) Copyright Nick Thompson 20182// (C) Copyright Matt Borland 20203// 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_DETAIL_SINGLE_PASS_HPP8#define BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_DETAIL_SINGLE_PASS_HPP9 10#include <boost/math/tools/config.hpp>11#include <boost/math/tools/assert.hpp>12#include <tuple>13#include <iterator>14#include <type_traits>15#include <cmath>16#include <algorithm>17#include <valarray>18#include <stdexcept>19#include <functional>20#include <vector>21 22#ifdef BOOST_MATH_HAS_THREADS23#include <future>24#include <thread>25#endif26 27namespace boost { namespace math { namespace statistics { namespace detail {28 29template<typename ReturnType, typename ForwardIterator>30ReturnType mean_sequential_impl(ForwardIterator first, ForwardIterator last)31{32 const std::size_t elements {static_cast<std::size_t>(std::distance(first, last))};33 std::valarray<ReturnType> mu {0, 0, 0, 0};34 std::valarray<ReturnType> temp {0, 0, 0, 0};35 ReturnType i {1};36 const ForwardIterator end {std::next(first, elements - (elements % 4))};37 ForwardIterator it {first};38 39 while(it != end)40 {41 const ReturnType inv {ReturnType(1) / i};42 temp = {static_cast<ReturnType>(*it++), static_cast<ReturnType>(*it++), static_cast<ReturnType>(*it++), static_cast<ReturnType>(*it++)};43 temp -= mu;44 mu += (temp *= inv);45 i += 1;46 }47 48 const ReturnType num1 {ReturnType(elements - (elements % 4))/ReturnType(4)};49 const ReturnType num2 {num1 + ReturnType(elements % 4)};50 51 while(it != last)52 {53 mu[3] += (*it-mu[3])/i;54 i += 1;55 ++it;56 }57 58 return (num1 * std::valarray<ReturnType>(mu[std::slice(0,3,1)]).sum() + num2 * mu[3]) / ReturnType(elements);59}60 61// Higham, Accuracy and Stability, equation 1.6a and 1.6b:62// Calculates Mean, M2, and variance63template<typename ReturnType, typename ForwardIterator>64ReturnType variance_sequential_impl(ForwardIterator first, ForwardIterator last)65{66 using Real = typename std::tuple_element<0, ReturnType>::type;67 68 Real M = *first;69 Real Q = 0;70 Real k = 2;71 Real M2 = 0;72 std::size_t n = 1;73 74 for(auto it = std::next(first); it != last; ++it)75 {76 Real tmp = (*it - M) / k;77 Real delta_1 = *it - M;78 Q += k*(k-1)*tmp*tmp;79 M += tmp;80 k += 1;81 Real delta_2 = *it - M;82 M2 += delta_1 * delta_2;83 ++n;84 }85 86 return std::make_tuple(M, M2, Q/(k-1), Real(n));87}88 89// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics90template<typename ReturnType, typename ForwardIterator>91ReturnType first_four_moments_sequential_impl(ForwardIterator first, ForwardIterator last)92{93 using Real = typename std::tuple_element<0, ReturnType>::type;94 using Size = typename std::tuple_element<4, ReturnType>::type;95 96 Real M1 = *first;97 Real M2 = 0;98 Real M3 = 0;99 Real M4 = 0;100 Size n = 2;101 for (auto it = std::next(first); it != last; ++it)102 {103 Real delta21 = *it - M1;104 Real tmp = delta21/n;105 M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);106 M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);107 M2 = M2 + tmp*(n-1)*delta21;108 M1 = M1 + tmp;109 n += 1;110 }111 112 return std::make_tuple(M1, M2, M3, M4, n-1);113}114 115#ifdef BOOST_MATH_HAS_THREADS116 117// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics118// EQN 3.1: https://www.osti.gov/servlets/purl/1426900119template<typename ReturnType, typename ForwardIterator>120ReturnType first_four_moments_parallel_impl(ForwardIterator first, ForwardIterator last)121{122 using Real = typename std::tuple_element<0, ReturnType>::type;123 124 const auto elements = std::distance(first, last);125 const unsigned max_concurrency = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();126 unsigned num_threads = 2u;127 128 // Threading is faster for: 10 + 5.13e-3 N/j <= 5.13e-3N => N >= 10^4j/5.13(j-1).129 const auto parallel_lower_bound = 10e4*max_concurrency/(5.13*(max_concurrency-1));130 const auto parallel_upper_bound = 10e4*2/5.13; // j = 2131 132 // https://lemire.me/blog/2020/01/30/cost-of-a-thread-in-c-under-linux/133 if(elements < parallel_lower_bound)134 {135 return detail::first_four_moments_sequential_impl<ReturnType>(first, last);136 }137 else if(elements >= parallel_upper_bound)138 {139 num_threads = max_concurrency;140 }141 else142 {143 for(unsigned i = 3; i < max_concurrency; ++i)144 {145 if(parallel_lower_bound < 10e4*i/(5.13*(i-1)))146 {147 num_threads = i;148 break;149 }150 }151 }152 153 std::vector<std::future<ReturnType>> future_manager;154 const auto elements_per_thread = std::ceil(static_cast<double>(elements) / num_threads);155 156 auto it = first;157 for(std::size_t i {}; i < num_threads - 1; ++i)158 {159 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, elements_per_thread]() -> ReturnType160 {161 return first_four_moments_sequential_impl<ReturnType>(it, std::next(it, elements_per_thread));162 }));163 it = std::next(it, elements_per_thread);164 }165 166 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, last]() -> ReturnType167 {168 return first_four_moments_sequential_impl<ReturnType>(it, last);169 }));170 171 auto temp = future_manager[0].get();172 Real M1_a = std::get<0>(temp);173 Real M2_a = std::get<1>(temp);174 Real M3_a = std::get<2>(temp);175 Real M4_a = std::get<3>(temp);176 Real range_a = std::get<4>(temp);177 178 for(std::size_t i = 1; i < future_manager.size(); ++i)179 {180 temp = future_manager[i].get();181 Real M1_b = std::get<0>(temp);182 Real M2_b = std::get<1>(temp);183 Real M3_b = std::get<2>(temp);184 Real M4_b = std::get<3>(temp);185 Real range_b = std::get<4>(temp);186 187 const Real n_ab = range_a + range_b;188 const Real delta = M1_b - M1_a;189 190 M1_a = (range_a * M1_a + range_b * M1_b) / n_ab;191 M2_a = M2_a + M2_b + delta * delta * (range_a * range_b / n_ab);192 M3_a = M3_a + M3_b + (delta * delta * delta) * range_a * range_b * (range_a - range_b) / (n_ab * n_ab) 193 + Real(3) * delta * (range_a * M2_b - range_b * M2_a) / n_ab;194 M4_a = M4_a + M4_b + (delta * delta * delta * delta) * range_a * range_b * (range_a * range_a - range_a * range_b + range_b * range_b) / (n_ab * n_ab * n_ab)195 + Real(6) * delta * delta * (range_a * range_a * M2_b + range_b * range_b * M2_a) / (n_ab * n_ab) 196 + Real(4) * delta * (range_a * M3_b - range_b * M3_a) / n_ab;197 range_a = n_ab;198 }199 200 return std::make_tuple(M1_a, M2_a, M3_a, M4_a, elements);201}202 203#endif // BOOST_MATH_HAS_THREADS204 205// Follows equation 1.5 of:206// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf207template<typename ReturnType, typename ForwardIterator>208ReturnType skewness_sequential_impl(ForwardIterator first, ForwardIterator last)209{210 using std::sqrt;211 BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute skewness.");212 213 ReturnType M1 = *first;214 ReturnType M2 = 0;215 ReturnType M3 = 0;216 ReturnType n = 2;217 218 for (auto it = std::next(first); it != last; ++it) 219 {220 ReturnType delta21 = *it - M1;221 ReturnType tmp = delta21/n;222 M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);223 M2 += tmp*(n-1)*delta21;224 M1 += tmp;225 n += 1;226 }227 228 ReturnType var = M2/(n-1);229 230 if (var == 0)231 {232 // The limit is technically undefined, but the interpretation here is clear:233 // A constant dataset has no skewness.234 return ReturnType(0);235 }236 237 ReturnType skew = M3/(M2*sqrt(var));238 return skew;239}240 241template<typename ReturnType, typename ForwardIterator>242ReturnType gini_coefficient_sequential_impl(ForwardIterator first, ForwardIterator last)243{244 ReturnType i = 1;245 ReturnType num = 0;246 ReturnType denom = 0;247 248 for(auto it = first; it != last; ++it)249 {250 num += *it*i;251 denom += *it;252 ++i;253 }254 255 // If the l1 norm is zero, all elements are zero, so every element is the same.256 if(denom == 0)257 {258 return ReturnType(0);259 }260 else261 {262 return ((2*num)/denom - i)/(i-1);263 }264}265 266template<typename ReturnType, typename ForwardIterator>267ReturnType gini_range_fraction(ForwardIterator first, ForwardIterator last, std::size_t starting_index)268{269 using Real = typename std::tuple_element<0, ReturnType>::type;270 271 std::size_t i = starting_index + 1;272 Real num = 0;273 Real denom = 0;274 275 for(auto it = first; it != last; ++it)276 {277 num += *it*i;278 denom += *it;279 ++i;280 }281 282 return std::make_tuple(num, denom, i);283}284 285#ifdef BOOST_MATH_HAS_THREADS286 287template<typename ReturnType, typename ExecutionPolicy, typename ForwardIterator>288ReturnType gini_coefficient_parallel_impl(ExecutionPolicy&&, ForwardIterator first, ForwardIterator last)289{290 using range_tuple = std::tuple<ReturnType, ReturnType, std::size_t>;291 292 const auto elements = std::distance(first, last);293 const unsigned max_concurrency = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();294 unsigned num_threads = 2u;295 296 // Threading is faster for: 10 + 10.12e-3 N/j <= 10.12e-3N => N >= 10^4j/10.12(j-1).297 const auto parallel_lower_bound = 10e4*max_concurrency/(10.12*(max_concurrency-1));298 const auto parallel_upper_bound = 10e4*2/10.12; // j = 2299 300 // https://lemire.me/blog/2020/01/30/cost-of-a-thread-in-c-under-linux/301 if(elements < parallel_lower_bound)302 {303 return gini_coefficient_sequential_impl<ReturnType>(first, last);304 }305 else if(elements >= parallel_upper_bound)306 {307 num_threads = max_concurrency;308 }309 else310 {311 for(unsigned i = 3; i < max_concurrency; ++i)312 {313 if(parallel_lower_bound < 10e4*i/(10.12*(i-1)))314 {315 num_threads = i;316 break;317 }318 }319 }320 321 std::vector<std::future<range_tuple>> future_manager;322 const auto elements_per_thread = std::ceil(static_cast<double>(elements) / num_threads);323 324 auto it = first;325 for(std::size_t i {}; i < num_threads - 1; ++i)326 {327 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, elements_per_thread, i]() -> range_tuple328 {329 return gini_range_fraction<range_tuple>(it, std::next(it, elements_per_thread), i*elements_per_thread);330 }));331 it = std::next(it, elements_per_thread);332 }333 334 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [it, last, num_threads, elements_per_thread]() -> range_tuple335 {336 return gini_range_fraction<range_tuple>(it, last, (num_threads - 1)*elements_per_thread);337 }));338 339 ReturnType num = 0;340 ReturnType denom = 0;341 342 for(std::size_t i = 0; i < future_manager.size(); ++i)343 {344 auto temp = future_manager[i].get();345 num += std::get<0>(temp);346 denom += std::get<1>(temp);347 }348 349 // If the l1 norm is zero, all elements are zero, so every element is the same.350 if(denom == 0)351 {352 return ReturnType(0);353 }354 else355 {356 return ((2*num)/denom - elements)/(elements-1);357 }358}359 360#endif // BOOST_MATH_HAS_THREADS361 362template<typename ForwardIterator, typename OutputIterator>363OutputIterator mode_impl(ForwardIterator first, ForwardIterator last, OutputIterator output)364{365 using Z = typename std::iterator_traits<ForwardIterator>::value_type;366 using Size = typename std::iterator_traits<ForwardIterator>::difference_type;367 368 std::vector<Z> modes {};369 modes.reserve(16);370 Size max_counter {0};371 372 while(first != last)373 {374 Size current_count {0};375 ForwardIterator end_it {first};376 while(end_it != last && *end_it == *first)377 {378 ++current_count;379 ++end_it;380 }381 382 if(current_count > max_counter)383 {384 modes.resize(1);385 modes[0] = *first;386 max_counter = current_count;387 }388 389 else if(current_count == max_counter)390 {391 modes.emplace_back(*first);392 }393 394 first = end_it;395 }396 397 return std::move(modes.begin(), modes.end(), output);398}399}}}}400 401#endif // BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_DETAIL_SINGLE_PASS_HPP402