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>
test_covariance()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>
test_correlation_coefficient()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
main()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 }
171