471 lines · plain
1// (C) Copyright Nick Thompson 2018.2// (C) Copyright Matt Borland 2021.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_BIVARIATE_STATISTICS_HPP8#define BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP9 10#include <iterator>11#include <tuple>12#include <type_traits>13#include <stdexcept>14#include <vector>15#include <algorithm>16#include <cmath>17#include <cstddef>18#include <boost/math/tools/assert.hpp>19#include <boost/math/tools/config.hpp>20 21#ifdef BOOST_MATH_EXEC_COMPATIBLE22#include <execution>23#include <future>24#include <thread>25#endif26 27namespace boost{ namespace math{ namespace statistics { namespace detail {28 29// See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.30template<typename ReturnType, typename ForwardIterator>31ReturnType means_and_covariance_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)32{33 using Real = typename std::tuple_element<0, ReturnType>::type;34 35 Real cov = 0;36 ForwardIterator u_it = u_begin;37 ForwardIterator v_it = v_begin;38 Real mu_u = *u_it++;39 Real mu_v = *v_it++;40 std::size_t i = 1;41 42 while(u_it != u_end && v_it != v_end)43 {44 Real u_temp = (*u_it++ - mu_u)/(i+1);45 Real v_temp = *v_it++ - mu_v;46 cov += i*u_temp*v_temp;47 mu_u = mu_u + u_temp;48 mu_v = mu_v + v_temp/(i+1);49 i = i + 1;50 }51 52 if(u_it != u_end || v_it != v_end)53 {54 throw std::domain_error("The size of each sample set must be the same to compute covariance");55 }56 57 return std::make_tuple(mu_u, mu_v, cov/i, Real(i));58}59 60#ifdef BOOST_MATH_EXEC_COMPATIBLE61 62// Numerically stable parallel computation of (co-)variance63// https://dl.acm.org/doi/10.1145/3221269.322303664template<typename ReturnType, typename ForwardIterator>65ReturnType means_and_covariance_parallel_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)66{67 using Real = typename std::tuple_element<0, ReturnType>::type;68 69 const auto u_elements = std::distance(u_begin, u_end);70 const auto v_elements = std::distance(v_begin, v_end);71 72 if(u_elements != v_elements)73 {74 throw std::domain_error("The size of each sample set must be the same to compute covariance");75 }76 77 const unsigned max_concurrency = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();78 unsigned num_threads = 2u;79 80 // 5.16 comes from benchmarking. See boost/math/reporting/performance/bivariate_statistics_performance.cpp81 // Threading is faster for: 10 + 5.16e-3 N/j <= 5.16e-3N => N >= 10^4j/5.16(j-1).82 const auto parallel_lower_bound = 10e4*max_concurrency/(5.16*(max_concurrency-1));83 const auto parallel_upper_bound = 10e4*2/5.16; // j = 284 85 // https://lemire.me/blog/2020/01/30/cost-of-a-thread-in-c-under-linux/86 if(u_elements < parallel_lower_bound)87 {88 return means_and_covariance_seq_impl<ReturnType>(u_begin, u_end, v_begin, v_end);89 }90 else if(u_elements >= parallel_upper_bound)91 {92 num_threads = max_concurrency;93 }94 else95 {96 for(unsigned i = 3; i < max_concurrency; ++i)97 {98 if(parallel_lower_bound < 10e4*i/(5.16*(i-1)))99 {100 num_threads = i;101 break;102 }103 }104 }105 106 std::vector<std::future<ReturnType>> future_manager;107 const auto elements_per_thread = std::ceil(static_cast<double>(u_elements)/num_threads);108 109 ForwardIterator u_it = u_begin;110 ForwardIterator v_it = v_begin;111 112 for(std::size_t i = 0; i < num_threads - 1; ++i)113 {114 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, v_it, elements_per_thread]() -> ReturnType115 {116 return means_and_covariance_seq_impl<ReturnType>(u_it, std::next(u_it, elements_per_thread), v_it, std::next(v_it, elements_per_thread));117 }));118 u_it = std::next(u_it, elements_per_thread);119 v_it = std::next(v_it, elements_per_thread);120 }121 122 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, u_end, v_it, v_end]() -> ReturnType123 {124 return means_and_covariance_seq_impl<ReturnType>(u_it, u_end, v_it, v_end);125 }));126 127 ReturnType temp = future_manager[0].get();128 Real mu_u_a = std::get<0>(temp);129 Real mu_v_a = std::get<1>(temp);130 Real cov_a = std::get<2>(temp);131 Real n_a = std::get<3>(temp);132 133 for(std::size_t i = 1; i < future_manager.size(); ++i)134 {135 temp = future_manager[i].get();136 Real mu_u_b = std::get<0>(temp);137 Real mu_v_b = std::get<1>(temp);138 Real cov_b = std::get<2>(temp);139 Real n_b = std::get<3>(temp);140 141 const Real n_ab = n_a + n_b;142 const Real delta_u = mu_u_b - mu_u_a;143 const Real delta_v = mu_v_b - mu_v_a;144 145 cov_a = cov_a + cov_b + (-delta_u)*(-delta_v)*((n_a*n_b)/n_ab);146 mu_u_a = mu_u_a + delta_u*(n_b/n_ab);147 mu_v_a = mu_v_a + delta_v*(n_b/n_ab);148 n_a = n_ab;149 }150 151 return std::make_tuple(mu_u_a, mu_v_a, cov_a, n_a);152}153 154#endif // BOOST_MATH_EXEC_COMPATIBLE155 156template<typename ReturnType, typename ForwardIterator>157ReturnType correlation_coefficient_seq_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)158{159 using Real = typename std::tuple_element<0, ReturnType>::type;160 using std::sqrt;161 162 Real cov = 0;163 ForwardIterator u_it = u_begin;164 ForwardIterator v_it = v_begin;165 Real mu_u = *u_it++;166 Real mu_v = *v_it++;167 Real Qu = 0;168 Real Qv = 0;169 std::size_t i = 1;170 171 while(u_it != u_end && v_it != v_end)172 {173 Real u_tmp = *u_it++ - mu_u;174 Real v_tmp = *v_it++ - mu_v;175 Qu = Qu + (i*u_tmp*u_tmp)/(i+1);176 Qv = Qv + (i*v_tmp*v_tmp)/(i+1);177 cov += i*u_tmp*v_tmp/(i+1);178 mu_u = mu_u + u_tmp/(i+1);179 mu_v = mu_v + v_tmp/(i+1);180 ++i;181 }182 183 184 // If one dataset is constant, then the correlation coefficient is undefined.185 // See https://stats.stackexchange.com/questions/23676/normalized-correlation-with-a-constant-vector186 // Thanks to zbjornson for pointing this out.187 if (Qu == 0 || Qv == 0)188 {189 return std::make_tuple(mu_u, Qu, mu_v, Qv, cov, std::numeric_limits<Real>::quiet_NaN(), Real(i));190 }191 192 // Make sure rho in [-1, 1], even in the presence of numerical noise.193 Real rho = cov/sqrt(Qu*Qv);194 if (rho > 1) {195 rho = 1;196 }197 if (rho < -1) {198 rho = -1;199 }200 201 return std::make_tuple(mu_u, Qu, mu_v, Qv, cov, rho, Real(i));202}203 204#ifdef BOOST_MATH_EXEC_COMPATIBLE205 206// Numerically stable parallel computation of (co-)variance:207// https://dl.acm.org/doi/10.1145/3221269.3223036208//209// Parallel computation of variance:210// http://i.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf211template<typename ReturnType, typename ForwardIterator>212ReturnType correlation_coefficient_parallel_impl(ForwardIterator u_begin, ForwardIterator u_end, ForwardIterator v_begin, ForwardIterator v_end)213{214 using Real = typename std::tuple_element<0, ReturnType>::type;215 216 const auto u_elements = std::distance(u_begin, u_end);217 const auto v_elements = std::distance(v_begin, v_end);218 219 if(u_elements != v_elements)220 {221 throw std::domain_error("The size of each sample set must be the same to compute covariance");222 }223 224 const unsigned max_concurrency = std::thread::hardware_concurrency() == 0 ? 2u : std::thread::hardware_concurrency();225 unsigned num_threads = 2u;226 227 // 3.25 comes from benchmarking. See boost/math/reporting/performance/bivariate_statistics_performance.cpp228 // Threading is faster for: 10 + 3.25e-3 N/j <= 3.25e-3N => N >= 10^4j/3.25(j-1).229 const auto parallel_lower_bound = 10e4*max_concurrency/(3.25*(max_concurrency-1));230 const auto parallel_upper_bound = 10e4*2/3.25; // j = 2231 232 // https://lemire.me/blog/2020/01/30/cost-of-a-thread-in-c-under-linux/233 if(u_elements < parallel_lower_bound)234 {235 return correlation_coefficient_seq_impl<ReturnType>(u_begin, u_end, v_begin, v_end);236 }237 else if(u_elements >= parallel_upper_bound)238 {239 num_threads = max_concurrency;240 }241 else242 {243 for(unsigned i = 3; i < max_concurrency; ++i)244 {245 if(parallel_lower_bound < 10e4*i/(3.25*(i-1)))246 {247 num_threads = i;248 break;249 }250 }251 }252 253 std::vector<std::future<ReturnType>> future_manager;254 const auto elements_per_thread = std::ceil(static_cast<double>(u_elements)/num_threads);255 256 ForwardIterator u_it = u_begin;257 ForwardIterator v_it = v_begin;258 259 for(std::size_t i = 0; i < num_threads - 1; ++i)260 {261 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, v_it, elements_per_thread]() -> ReturnType262 {263 return correlation_coefficient_seq_impl<ReturnType>(u_it, std::next(u_it, elements_per_thread), v_it, std::next(v_it, elements_per_thread));264 }));265 u_it = std::next(u_it, elements_per_thread);266 v_it = std::next(v_it, elements_per_thread);267 }268 269 future_manager.emplace_back(std::async(std::launch::async | std::launch::deferred, [u_it, u_end, v_it, v_end]() -> ReturnType270 {271 return correlation_coefficient_seq_impl<ReturnType>(u_it, u_end, v_it, v_end);272 }));273 274 ReturnType temp = future_manager[0].get();275 Real mu_u_a = std::get<0>(temp);276 Real Qu_a = std::get<1>(temp);277 Real mu_v_a = std::get<2>(temp);278 Real Qv_a = std::get<3>(temp);279 Real cov_a = std::get<4>(temp);280 Real n_a = std::get<6>(temp);281 282 for(std::size_t i = 1; i < future_manager.size(); ++i)283 {284 temp = future_manager[i].get();285 Real mu_u_b = std::get<0>(temp);286 Real Qu_b = std::get<1>(temp);287 Real mu_v_b = std::get<2>(temp);288 Real Qv_b = std::get<3>(temp);289 Real cov_b = std::get<4>(temp);290 Real n_b = std::get<6>(temp);291 292 const Real n_ab = n_a + n_b;293 const Real delta_u = mu_u_b - mu_u_a;294 const Real delta_v = mu_v_b - mu_v_a;295 296 cov_a = cov_a + cov_b + (-delta_u)*(-delta_v)*((n_a*n_b)/n_ab);297 mu_u_a = mu_u_a + delta_u*(n_b/n_ab);298 mu_v_a = mu_v_a + delta_v*(n_b/n_ab);299 Qu_a = Qu_a + Qu_b + delta_u*delta_u*((n_a*n_b)/n_ab);300 Qv_b = Qv_a + Qv_b + delta_v*delta_v*((n_a*n_b)/n_ab);301 n_a = n_ab;302 }303 304 // If one dataset is constant, then the correlation coefficient is undefined.305 // See https://stats.stackexchange.com/questions/23676/normalized-correlation-with-a-constant-vector306 // Thanks to zbjornson for pointing this out.307 if (Qu_a == 0 || Qv_a == 0)308 {309 return std::make_tuple(mu_u_a, Qu_a, mu_v_a, Qv_a, cov_a, std::numeric_limits<Real>::quiet_NaN(), n_a);310 }311 312 // Make sure rho in [-1, 1], even in the presence of numerical noise.313 Real rho = cov_a/sqrt(Qu_a*Qv_a);314 if (rho > 1) {315 rho = 1;316 }317 if (rho < -1) {318 rho = -1;319 }320 321 return std::make_tuple(mu_u_a, Qu_a, mu_v_a, Qv_a, cov_a, rho, n_a);322}323 324#endif // BOOST_MATH_EXEC_COMPATIBLE325 326} // namespace detail327 328#ifdef BOOST_MATH_EXEC_COMPATIBLE329 330template<typename ExecutionPolicy, typename Container, typename Real = typename Container::value_type>331inline auto means_and_covariance(ExecutionPolicy&& exec, Container const & u, Container const & v)332{333 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)334 {335 if constexpr (std::is_integral_v<Real>)336 {337 using ReturnType = std::tuple<double, double, double, double>;338 ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));339 return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));340 }341 else342 {343 using ReturnType = std::tuple<Real, Real, Real, Real>;344 ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));345 return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));346 }347 }348 else349 {350 if constexpr (std::is_integral_v<Real>)351 {352 using ReturnType = std::tuple<double, double, double, double>;353 ReturnType temp = detail::means_and_covariance_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));354 return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));355 }356 else357 {358 using ReturnType = std::tuple<Real, Real, Real, Real>;359 ReturnType temp = detail::means_and_covariance_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));360 return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));361 }362 }363}364 365template<typename Container>366inline auto means_and_covariance(Container const & u, Container const & v)367{368 return means_and_covariance(std::execution::seq, u, v);369}370 371template<typename ExecutionPolicy, typename Container>372inline auto covariance(ExecutionPolicy&& exec, Container const & u, Container const & v)373{374 return std::get<2>(means_and_covariance(exec, u, v));375}376 377template<typename Container>378inline auto covariance(Container const & u, Container const & v)379{380 return covariance(std::execution::seq, u, v);381}382 383template<typename ExecutionPolicy, typename Container, typename Real = typename Container::value_type>384inline auto correlation_coefficient(ExecutionPolicy&& exec, Container const & u, Container const & v)385{386 if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>)387 {388 if constexpr (std::is_integral_v<Real>)389 {390 using ReturnType = std::tuple<double, double, double, double, double, double, double>;391 return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));392 }393 else394 {395 using ReturnType = std::tuple<Real, Real, Real, Real, Real, Real, Real>;396 return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));397 }398 }399 else400 {401 if constexpr (std::is_integral_v<Real>)402 {403 using ReturnType = std::tuple<double, double, double, double, double, double, double>;404 return std::get<5>(detail::correlation_coefficient_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));405 }406 else407 {408 using ReturnType = std::tuple<Real, Real, Real, Real, Real, Real, Real>;409 return std::get<5>(detail::correlation_coefficient_parallel_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));410 }411 }412}413 414template<typename Container, typename Real = typename Container::value_type>415inline auto correlation_coefficient(Container const & u, Container const & v)416{417 return correlation_coefficient(std::execution::seq, u, v);418}419 420#else // C++11 and single threaded bindings421 422template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>423inline auto means_and_covariance(Container const & u, Container const & v) -> std::tuple<double, double, double>424{425 using ReturnType = std::tuple<double, double, double, double>;426 ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));427 return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));428}429 430template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>431inline auto means_and_covariance(Container const & u, Container const & v) -> std::tuple<Real, Real, Real>432{433 using ReturnType = std::tuple<Real, Real, Real, Real>;434 ReturnType temp = detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v));435 return std::make_tuple(std::get<0>(temp), std::get<1>(temp), std::get<2>(temp));436}437 438template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>439inline double covariance(Container const & u, Container const & v)440{441 using ReturnType = std::tuple<double, double, double, double>;442 return std::get<2>(detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));443}444 445template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>446inline Real covariance(Container const & u, Container const & v)447{448 using ReturnType = std::tuple<Real, Real, Real, Real>;449 return std::get<2>(detail::means_and_covariance_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));450}451 452template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>453inline double correlation_coefficient(Container const & u, Container const & v)454{455 using ReturnType = std::tuple<double, double, double, double, double, double, double>;456 return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));457}458 459template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>460inline Real correlation_coefficient(Container const & u, Container const & v)461{462 using ReturnType = std::tuple<Real, Real, Real, Real, Real, Real, Real>;463 return std::get<5>(detail::correlation_coefficient_seq_impl<ReturnType>(std::begin(u), std::end(u), std::begin(v), std::end(v)));464}465 466#endif467 468}}} // namespace boost::math::statistics469 470#endif471