boost/math/tools/bivariate_statistics.hpp

   1 //  (C) Copyright Nick Thompson 2018.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_BIVARIATE_STATISTICS_HPP
   7 #define BOOST_MATH_TOOLS_BIVARIATE_STATISTICS_HPP
   8
   9 #include <iterator>
  10 #include <tuple>
  11 #include <boost/assert.hpp>
  12 #include <boost/config/header_deprecated.hpp>
  13
  14 BOOST_HEADER_DEPRECATED("<boost/math/statistics/bivariate_statistics.hpp>");
  15
  16 namespace boost{ namespace math{ namespace tools {
  17
  18 template<class Container>
  19 auto means_and_covariance(Container const & u, Container const & v)
  20 {
  21     using Real = typename Container::value_type;
  22     using std::size;
  23     BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
  24     BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least one sample.");
  25
  26     // See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.
  27     Real cov = 0;
  28     Real mu_u = u[0];
  29     Real mu_v = v[0];
  30
  31     for(size_t i = 1; i < size(u); ++i)
  32     {
  33         Real u_tmp = (u[i] - mu_u)/(i+1);
  34         Real v_tmp = v[i] - mu_v;
  35         cov += i*u_tmp*v_tmp;
  36         mu_u = mu_u + u_tmp;
  37         mu_v = mu_v + v_tmp/(i+1);
  38     }
  39
  40     return std::make_tuple(mu_u, mu_v, cov/size(u));
  41 }
  42
  43 template<class Container>
  44 auto covariance(Container const & u, Container const & v)
  45 {
  46     auto [mu_u, mu_v, cov] = boost::math::tools::means_and_covariance(u, v);
  47     return cov;
  48 }
  49
  50 template<class Container>
  51 auto correlation_coefficient(Container const & u, Container const & v)
  52 {
  53     using Real = typename Container::value_type;
  54     using std::size;
  55     BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
  56     BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least two samples.");
  57
  58     Real cov = 0;
  59     Real mu_u = u[0];
  60     Real mu_v = v[0];
  61     Real Qu = 0;
  62     Real Qv = 0;
  63
  64     for(size_t i = 1; i < size(u); ++i)
  65     {
  66         Real u_tmp = u[i] - mu_u;
  67         Real v_tmp = v[i] - mu_v;
  68         Qu = Qu + (i*u_tmp*u_tmp)/(i+1);
  69         Qv = Qv + (i*v_tmp*v_tmp)/(i+1);
  70         cov += i*u_tmp*v_tmp/(i+1);
  71         mu_u = mu_u + u_tmp/(i+1);
  72         mu_v = mu_v + v_tmp/(i+1);
  73     }
  74
  75     // If both datasets are constant, then they are perfectly correlated.
  76     if (Qu == 0 && Qv == 0)
  77     {
  78         return Real(1);
  79     }
  80     // If one dataset is constant and the other isn't, then they have no correlation:
  81     if (Qu == 0 || Qv == 0)
  82     {
  83         return Real(0);
  84     }
  85
  86     // Make sure rho in [-1, 1], even in the presence of numerical noise.
  87     Real rho = cov/sqrt(Qu*Qv);
  88     if (rho > 1) {
  89         rho = 1;
  90     }
  91     if (rho < -1) {
  92         rho = -1;
  93     }
  94     return rho;
  95 }
  96
  97 }}}
  98 #endif