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)
6 #ifndef BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP
7 #define BOOST_MATH_STATISTICS_BIVARIATE_STATISTICS_HPP
11 #include <boost/assert.hpp>
14 namespace boost{ namespace math{ namespace statistics {
16 template<class Container>
17 auto means_and_covariance(Container const & u, Container const & v)
19 using Real = typename Container::value_type;
21 BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
22 BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least one sample.");
24 // See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.
29 for(size_t i = 1; i < size(u); ++i)
31 Real u_tmp = (u[i] - mu_u)/(i+1);
32 Real v_tmp = v[i] - mu_v;
35 mu_v = mu_v + v_tmp/(i+1);
38 return std::make_tuple(mu_u, mu_v, cov/size(u));
41 template<class Container>
42 auto covariance(Container const & u, Container const & v)
44 auto [mu_u, mu_v, cov] = boost::math::statistics::means_and_covariance(u, v);
48 template<class Container>
49 auto correlation_coefficient(Container const & u, Container const & v)
51 using Real = typename Container::value_type;
53 BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
54 BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least two samples.");
62 for(size_t i = 1; i < size(u); ++i)
64 Real u_tmp = u[i] - mu_u;
65 Real v_tmp = v[i] - mu_v;
66 Qu = Qu + (i*u_tmp*u_tmp)/(i+1);
67 Qv = Qv + (i*v_tmp*v_tmp)/(i+1);
68 cov += i*u_tmp*v_tmp/(i+1);
69 mu_u = mu_u + u_tmp/(i+1);
70 mu_v = mu_v + v_tmp/(i+1);
73 // If both datasets are constant, then they are perfectly correlated.
74 if (Qu == 0 && Qv == 0)
78 // If one dataset is constant and the other isn't, then they have no correlation:
79 if (Qu == 0 || Qv == 0)
84 // Make sure rho in [-1, 1], even in the presence of numerical noise.
85 Real rho = cov/sqrt(Qu*Qv);