97 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_BIVARIATE_STATISTICS_HPP7#define BOOST_MATH_TOOLS_BIVARIATE_STATISTICS_HPP8 9#include <iterator>10#include <tuple>11#include <limits>12#include <boost/math/tools/assert.hpp>13#include <boost/math/tools/header_deprecated.hpp>14 15BOOST_MATH_HEADER_DEPRECATED("<boost/math/statistics/bivariate_statistics.hpp>");16 17namespace boost{ namespace math{ namespace tools {18 19template<class Container>20auto means_and_covariance(Container const & u, Container const & v)21{22 using Real = typename Container::value_type;23 using std::size;24 BOOST_MATH_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");25 BOOST_MATH_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least one sample.");26 27 // See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.28 Real cov = 0;29 Real mu_u = u[0];30 Real mu_v = v[0];31 32 for(size_t i = 1; i < size(u); ++i)33 {34 Real u_tmp = (u[i] - mu_u)/(i+1);35 Real v_tmp = v[i] - mu_v;36 cov += i*u_tmp*v_tmp;37 mu_u = mu_u + u_tmp;38 mu_v = mu_v + v_tmp/(i+1);39 }40 41 return std::make_tuple(mu_u, mu_v, cov/size(u));42}43 44template<class Container>45auto covariance(Container const & u, Container const & v)46{47 auto [mu_u, mu_v, cov] = boost::math::tools::means_and_covariance(u, v);48 return cov;49}50 51template<class Container>52auto correlation_coefficient(Container const & u, Container const & v)53{54 using Real = typename Container::value_type;55 using std::size;56 BOOST_MATH_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");57 BOOST_MATH_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least two samples.");58 59 Real cov = 0;60 Real mu_u = u[0];61 Real mu_v = v[0];62 Real Qu = 0;63 Real Qv = 0;64 65 for(size_t i = 1; i < size(u); ++i)66 {67 Real u_tmp = u[i] - mu_u;68 Real v_tmp = v[i] - mu_v;69 Qu = Qu + (i*u_tmp*u_tmp)/(i+1);70 Qv = Qv + (i*v_tmp*v_tmp)/(i+1);71 cov += i*u_tmp*v_tmp/(i+1);72 mu_u = mu_u + u_tmp/(i+1);73 mu_v = mu_v + v_tmp/(i+1);74 }75 76 // If one dataset is constant, then they have no correlation:77 // See https://stats.stackexchange.com/questions/23676/normalized-correlation-with-a-constant-vector78 // Thanks to zbjornson for pointing this out.79 if (Qu == 0 || Qv == 0)80 {81 return std::numeric_limits<Real>::quiet_NaN();82 }83 84 // Make sure rho in [-1, 1], even in the presence of numerical noise.85 Real rho = cov/sqrt(Qu*Qv);86 if (rho > 1) {87 rho = 1;88 }89 if (rho < -1) {90 rho = -1;91 }92 return rho;93}94 95}}}96#endif97