libs/math/test/bivariate_statistics_test.cpp

   1 /*
   2  *  (C) Copyright Nick Thompson 2018.
   3  *  Use, modification and distribution are subject to the
   4  *  Boost Software License, Version 1.0. (See accompanying file
   5  *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #include <iostream>
   9 #include <iomanip>
  10 #include <vector>
  11 #include <array>
  12 #include <forward_list>
  13 #include <algorithm>
  14 #include <random>
  15 #include <boost/core/lightweight_test.hpp>
  16 #include <boost/numeric/ublas/vector.hpp>
  17 #include <boost/math/constants/constants.hpp>
  18 #include <boost/math/statistics/univariate_statistics.hpp>
  19 #include <boost/math/statistics/bivariate_statistics.hpp>
  20 #include <boost/multiprecision/cpp_bin_float.hpp>
  21 #include <boost/multiprecision/cpp_complex.hpp>
  22
  23 using boost::multiprecision::cpp_bin_float_50;
  24 using boost::multiprecision::cpp_complex_50;
  25
  26 /*
  27  * Test checklist:
  28  * 1) Does it work with multiprecision?
  29  * 2) Does it work with .cbegin()/.cend() if the data is not altered?
  30  * 3) Does it work with ublas and std::array? (Checking Eigen and Armadillo will make the CI system really unhappy.)
  31  * 4) Does it work with std::forward_list if a forward iterator is all that is required?
  32  * 5) Does it work with complex data if complex data is sensible?
  33  */
  34
  35 using  boost::math::statistics::means_and_covariance;
  36 using  boost::math::statistics::covariance;
  37
  38 template<class Real>
  39 void test_covariance()
  40 {
  41     std::cout << std::setprecision(std::numeric_limits<Real>::digits10+1);
  42     Real tol = std::numeric_limits<Real>::epsilon();
  43     using std::abs;
  44
  45     // Covariance of a single thing is zero:
  46     std::array<Real, 1> u1{8};
  47     std::array<Real, 1> v1{17};
  48     auto [mu_u1, mu_v1, cov1] = means_and_covariance(u1, v1);
  49
  50     BOOST_TEST(abs(cov1) < tol);
  51     BOOST_TEST(abs(mu_u1 - 8) < tol);
  52     BOOST_TEST(abs(mu_v1 - 17) < tol);
  53
  54
  55     std::array<Real, 2> u2{8, 4};
  56     std::array<Real, 2> v2{3, 7};
  57     auto [mu_u2, mu_v2, cov2] = means_and_covariance(u2, v2);
  58
  59     BOOST_TEST(abs(cov2+4) < tol);
  60     BOOST_TEST(abs(mu_u2 - 6) < tol);
  61     BOOST_TEST(abs(mu_v2 - 5) < tol);
  62
  63     std::vector<Real> u3{1,2,3};
  64     std::vector<Real> v3{1,1,1};
  65
  66     auto [mu_u3, mu_v3, cov3] = means_and_covariance(u3, v3);
  67
  68     // Since v is constant, covariance(u,v) = 0 against everything any u:
  69     BOOST_TEST(abs(cov3) < tol);
  70     BOOST_TEST(abs(mu_u3 - 2) < tol);
  71     BOOST_TEST(abs(mu_v3 - 1) < tol);
  72     // Make sure we pull the correct symbol out of means_and_covariance:
  73     cov3 = covariance(u3, v3);
  74     BOOST_TEST(abs(cov3) < tol);
  75
  76     cov3 = covariance(v3, u3);
  77     // Covariance is symmetric: cov(u,v) = cov(v,u)
  78     BOOST_TEST(abs(cov3) < tol);
  79
  80     // cov(u,u) = sigma(u)^2:
  81     cov3 = covariance(u3, u3);
  82     Real expected = Real(2)/Real(3);
  83
  84     BOOST_TEST(abs(cov3 - expected) < tol);
  85
  86     std::mt19937 gen(15);
  87     // Can't template standard library on multiprecision, so use double and cast back:
  88     std::uniform_real_distribution<double> dis(-1.0, 1.0);
  89     std::vector<Real> u(500);
  90     std::vector<Real> v(500);
  91     for(size_t i = 0; i < u.size(); ++i)
  92     {
  93         u[i] = (Real) dis(gen);
  94         v[i] = (Real) dis(gen);
  95     }
  96
  97     Real mu_u = boost::math::statistics::mean(u);
  98     Real mu_v = boost::math::statistics::mean(v);
  99     Real sigma_u_sq = boost::math::statistics::variance(u);
 100     Real sigma_v_sq = boost::math::statistics::variance(v);
 101
 102     auto [mu_u_, mu_v_, cov_uv] = means_and_covariance(u, v);
 103     BOOST_TEST(abs(mu_u - mu_u_) < tol);
 104     BOOST_TEST(abs(mu_v - mu_v_) < tol);
 105
 106     // Cauchy-Schwartz inequality:
 107     BOOST_TEST(cov_uv*cov_uv <= sigma_u_sq*sigma_v_sq);
 108     // cov(X, X) = sigma(X)^2:
 109     Real cov_uu = covariance(u, u);
 110     BOOST_TEST(abs(cov_uu - sigma_u_sq) < tol);
 111     Real cov_vv = covariance(v, v);
 112     BOOST_TEST(abs(cov_vv - sigma_v_sq) < tol);
 113
 114 }
 115
 116 template<class Real>
 117 void test_correlation_coefficient()
 118 {
 119     using boost::math::statistics::correlation_coefficient;
 120
 121     Real tol = std::numeric_limits<Real>::epsilon();
 122     std::vector<Real> u{1};
 123     std::vector<Real> v{1};
 124     Real rho_uv = correlation_coefficient(u, v);
 125     BOOST_TEST(abs(rho_uv - 1) < tol);
 126
 127     u = {1,1};
 128     v = {1,1};
 129     rho_uv = correlation_coefficient(u, v);
 130     BOOST_TEST(abs(rho_uv - 1) < tol);
 131
 132     u = {1, 2, 3};
 133     v = {1, 2, 3};
 134     rho_uv = correlation_coefficient(u, v);
 135     BOOST_TEST(abs(rho_uv - 1) < tol);
 136
 137     u = {1, 2, 3};
 138     v = {-1, -2, -3};
 139     rho_uv = correlation_coefficient(u, v);
 140     BOOST_TEST(abs(rho_uv + 1) < tol);
 141
 142     rho_uv = correlation_coefficient(v, u);
 143     BOOST_TEST(abs(rho_uv + 1) < tol);
 144
 145     u = {1, 2, 3};
 146     v = {0, 0, 0};
 147     rho_uv = correlation_coefficient(v, u);
 148     BOOST_TEST(abs(rho_uv) < tol);
 149
 150     u = {1, 2, 3};
 151     v = {0, 0, 3};
 152     rho_uv = correlation_coefficient(v, u);
 153     // mu_u = 2, sigma_u^2 = 2/3, mu_v = 1, sigma_v^2 = 2, cov(u,v) = 1.
 154     BOOST_TEST(abs(rho_uv - sqrt(Real(3))/Real(2)) < tol);
 155 }
 156
 157 int main()
 158 {
 159     test_covariance<float>();
 160     test_covariance<double>();
 161     test_covariance<long double>();
 162     test_covariance<cpp_bin_float_50>();
 163
 164     test_correlation_coefficient<float>();
 165     test_correlation_coefficient<double>();
 166     test_correlation_coefficient<long double>();
 167     test_correlation_coefficient<cpp_bin_float_50>();
 168
 169     return boost::report_errors();
 170 }