libs/math/test/signal_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 <vector>
   9 #include <array>
  10 #include <forward_list>
  11 #include <algorithm>
  12 #include <random>
  13 #include <boost/core/lightweight_test.hpp>
  14 #include <boost/numeric/ublas/vector.hpp>
  15 #include <boost/math/constants/constants.hpp>
  16 #include <boost/math/statistics/univariate_statistics.hpp>
  17 #include <boost/math/statistics/signal_statistics.hpp>
  18 #include <boost/multiprecision/cpp_bin_float.hpp>
  19 #include <boost/multiprecision/cpp_complex.hpp>
  20
  21 using std::abs;
  22 using boost::multiprecision::cpp_bin_float_50;
  23 using boost::multiprecision::cpp_complex_50;
  24 using boost::math::constants::two_pi;
  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  * 6) Does it work with integer data if sensible?
  34  */
  35
  36 template<class Real>
  37 void test_hoyer_sparsity()
  38 {
  39     using std::sqrt;
  40     Real tol = 5*std::numeric_limits<Real>::epsilon();
  41     std::vector<Real> v{1,0,0};
  42     Real hs = boost::math::statistics::hoyer_sparsity(v.begin(), v.end());
  43     BOOST_TEST(abs(hs - 1) < tol);
  44
  45     hs = boost::math::statistics::hoyer_sparsity(v);
  46     BOOST_TEST(abs(hs - 1) < tol);
  47
  48     // Does it work with constant iterators?
  49     hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend());
  50     BOOST_TEST(abs(hs - 1) < tol);
  51
  52     v[0] = 1;
  53     v[1] = 1;
  54     v[2] = 1;
  55     hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend());
  56     BOOST_TEST(abs(hs) < tol);
  57
  58     std::array<Real, 3> w{1,1,1};
  59     hs = boost::math::statistics::hoyer_sparsity(w);
  60     BOOST_TEST(abs(hs) < tol);
  61
  62     // Now some statistics:
  63     // If x_i ~ Unif(0,1), E[x_i] = 1/2, E[x_i^2] = 1/3.
  64     // Therefore, E[||x||_1] = N/2, E[||x||_2] = sqrt(N/3),
  65     // and hoyer_sparsity(x) is close to (1-sqrt(3)/2)/(1-1/sqrt(N))
  66     std::mt19937 gen(82);
  67     std::uniform_real_distribution<long double> dis(0, 1);
  68     v.resize(5000);
  69     for (size_t i = 0; i < v.size(); ++i) {
  70         v[i] = dis(gen);
  71     }
  72     hs = boost::math::statistics::hoyer_sparsity(v);
  73     Real expected = (1.0 - boost::math::constants::root_three<Real>()/2)/(1.0 - 1.0/sqrt(v.size()));
  74     BOOST_TEST(abs(expected - hs) < 0.01);
  75
  76     // Does it work with a forward list?
  77     std::forward_list<Real> u1{1, 1, 1};
  78     hs = boost::math::statistics::hoyer_sparsity(u1);
  79     BOOST_TEST(abs(hs) < tol);
  80
  81     // Does it work with a boost ublas vector?
  82     boost::numeric::ublas::vector<Real> u2(3);
  83     u2[0] = 1;
  84     u2[1] = 1;
  85     u2[2] = 1;
  86     hs = boost::math::statistics::hoyer_sparsity(u2);
  87     BOOST_TEST(abs(hs) < tol);
  88
  89 }
  90
  91 template<class Z>
  92 void test_integer_hoyer_sparsity()
  93 {
  94     using std::sqrt;
  95     double tol = 5*std::numeric_limits<double>::epsilon();
  96     std::vector<Z> v{1,0,0};
  97     double hs = boost::math::statistics::hoyer_sparsity(v);
  98     BOOST_TEST(abs(hs - 1) < tol);
  99
 100     v[0] = 1;
 101     v[1] = 1;
 102     v[2] = 1;
 103     hs = boost::math::statistics::hoyer_sparsity(v);
 104     BOOST_TEST(abs(hs) < tol);
 105 }
 106
 107
 108 template<class Complex>
 109 void test_complex_hoyer_sparsity()
 110 {
 111     typedef typename Complex::value_type Real;
 112     using std::sqrt;
 113     Real tol = 5*std::numeric_limits<Real>::epsilon();
 114     std::vector<Complex> v{{0,1}, {0, 0}, {0,0}};
 115     Real hs = boost::math::statistics::hoyer_sparsity(v.begin(), v.end());
 116     BOOST_TEST(abs(hs - 1) < tol);
 117
 118     hs = boost::math::statistics::hoyer_sparsity(v);
 119     BOOST_TEST(abs(hs - 1) < tol);
 120
 121     // Does it work with constant iterators?
 122     hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend());
 123     BOOST_TEST(abs(hs - 1) < tol);
 124
 125     // All are the same magnitude:
 126     v[0] = {0, 1};
 127     v[1] = {1, 0};
 128     v[2] = {0,-1};
 129     hs = boost::math::statistics::hoyer_sparsity(v.cbegin(), v.cend());
 130     BOOST_TEST(abs(hs) < tol);
 131 }
 132
 133
 134 template<class Real>
 135 void test_absolute_gini_coefficient()
 136 {
 137     using boost::math::statistics::absolute_gini_coefficient;
 138     using boost::math::statistics::sample_absolute_gini_coefficient;
 139     Real tol = std::numeric_limits<Real>::epsilon();
 140     std::vector<Real> v{-1,0,0};
 141     Real gini = sample_absolute_gini_coefficient(v.begin(), v.end());
 142     BOOST_TEST(abs(gini - 1) < tol);
 143
 144     gini = absolute_gini_coefficient(v);
 145     BOOST_TEST(abs(gini - Real(2)/Real(3)) < tol);
 146
 147     v[0] = 1;
 148     v[1] = -1;
 149     v[2] = 1;
 150     gini = absolute_gini_coefficient(v.begin(), v.end());
 151     BOOST_TEST(abs(gini) < tol);
 152     gini = sample_absolute_gini_coefficient(v.begin(), v.end());
 153     BOOST_TEST(abs(gini) < tol);
 154
 155     std::vector<std::complex<Real>> w(128);
 156     std::complex<Real> i{0,1};
 157     for(size_t k = 0; k < w.size(); ++k)
 158     {
 159         w[k] = exp(i*static_cast<Real>(k)/static_cast<Real>(w.size()));
 160     }
 161     gini = absolute_gini_coefficient(w.begin(), w.end());
 162     BOOST_TEST(abs(gini) < tol);
 163     gini = sample_absolute_gini_coefficient(w.begin(), w.end());
 164     BOOST_TEST(abs(gini) < tol);
 165
 166     // The population Gini index is invariant under "cloning": If w = v \oplus v, then G(w) = G(v).
 167     // We use the sample Gini index, so we need to rescale
 168     std::vector<Real> u(1000);
 169     std::mt19937 gen(35);
 170     std::uniform_real_distribution<long double> dis(0, 50);
 171     for (size_t i = 0; i < u.size()/2; ++i)
 172     {
 173         u[i] = dis(gen);
 174     }
 175     for (size_t i = 0; i < u.size()/2; ++i)
 176     {
 177         u[i + u.size()/2] = u[i];
 178     }
 179     Real population_gini1 = absolute_gini_coefficient(u.begin(), u.begin() + u.size()/2);
 180     Real population_gini2 = absolute_gini_coefficient(u.begin(), u.end());
 181
 182     BOOST_TEST(abs(population_gini1 - population_gini2) < 10*tol);
 183
 184     // The Gini coefficient of a uniform distribution is (b-a)/(3*(b+a)), see https://en.wikipedia.org/wiki/Gini_coefficient
 185     Real expected = (dis.b() - dis.a() )/(3*(dis.a() + dis.b()));
 186
 187     BOOST_TEST(abs(expected - population_gini1) < 0.01);
 188
 189     std::exponential_distribution<long double> exp_dis(1);
 190     for (size_t i = 0; i < u.size(); ++i)
 191     {
 192         u[i] = exp_dis(gen);
 193     }
 194     population_gini2 = absolute_gini_coefficient(u);
 195
 196     BOOST_TEST(abs(population_gini2 - 0.5) < 0.01);
 197 }
 198
 199
 200 template<class Real>
 201 void test_oracle_snr()
 202 {
 203     using std::abs;
 204     Real tol = 100*std::numeric_limits<Real>::epsilon();
 205     size_t length = 100;
 206     std::vector<Real> signal(length, 1);
 207     std::vector<Real> noisy_signal = signal;
 208
 209     noisy_signal[0] += 1;
 210     Real snr = boost::math::statistics::oracle_snr(signal, noisy_signal);
 211     Real snr_db = boost::math::statistics::oracle_snr_db(signal, noisy_signal);
 212     BOOST_TEST(abs(snr - length) < tol);
 213     BOOST_TEST(abs(snr_db - 10*log10(length)) < tol);
 214 }
 215
 216 template<class Z>
 217 void test_integer_oracle_snr()
 218 {
 219     double tol = std::numeric_limits<double>::epsilon();
 220     size_t length = 100;
 221     std::vector<Z> signal(length, 1);
 222     std::vector<Z> noisy_signal = signal;
 223
 224     noisy_signal[0] += 1;
 225     double snr = boost::math::statistics::oracle_snr(signal, noisy_signal);
 226     double snr_db = boost::math::statistics::oracle_snr_db(signal, noisy_signal);
 227     BOOST_TEST(abs(snr - length) < tol);
 228     BOOST_TEST(abs(snr_db - 10*log10(length)) < tol);
 229 }
 230
 231 template<class Complex>
 232 void test_complex_oracle_snr()
 233 {
 234     using Real = typename Complex::value_type;
 235     using std::abs;
 236     using std::log10;
 237     Real tol = 100*std::numeric_limits<Real>::epsilon();
 238     size_t length = 100;
 239     std::vector<Complex> signal(length, {1,0});
 240     std::vector<Complex> noisy_signal = signal;
 241
 242     noisy_signal[0] += Complex(1,0);
 243     Real snr = boost::math::statistics::oracle_snr(signal, noisy_signal);
 244     Real snr_db = boost::math::statistics::oracle_snr_db(signal, noisy_signal);
 245     BOOST_TEST(abs(snr - length) < tol);
 246     BOOST_TEST(abs(snr_db - 10*log10(length)) < tol);
 247 }
 248
 249 template<class Real>
 250 void test_m2m4_snr_estimator()
 251 {
 252     Real tol = std::numeric_limits<Real>::epsilon();
 253     std::vector<Real> signal(5000, 1);
 254     std::vector<Real> x(signal.size());
 255     std::mt19937 gen(18);
 256     std::normal_distribution<Real> dis{0, 1.0};
 257
 258     for (size_t i = 0; i < x.size(); ++i)
 259     {
 260         signal[i] = 5*sin(100*6.28*i/x.size());
 261         x[i] = signal[i] + dis(gen);
 262     }
 263
 264     // Kurtosis of a sine wave is 1.5:
 265     auto m2m4_db = boost::math::statistics::m2m4_snr_estimator_db(x, 1.5);
 266     auto oracle_snr_db = boost::math::statistics::mean_invariant_oracle_snr_db(signal, x);
 267     BOOST_TEST(abs(m2m4_db - oracle_snr_db) < 0.2);
 268
 269     std::uniform_real_distribution<Real> uni_dis{-1,1};
 270     for (size_t i = 0; i < x.size(); ++i)
 271     {
 272         x[i] = signal[i] + uni_dis(gen);
 273     }
 274
 275     // Kurtosis of continuous uniform distribution over [-1,1] is 1.8:
 276     m2m4_db = boost::math::statistics::m2m4_snr_estimator_db(x, 1.5, 1.8);
 277     oracle_snr_db = boost::math::statistics::mean_invariant_oracle_snr_db(signal, x);
 278     // The performance depends on the exact numbers generated by the distribution, but this isn't bad:
 279     BOOST_TEST(abs(m2m4_db - oracle_snr_db) < 0.2);
 280
 281     // The SNR estimator should be scale invariant.
 282     // If x has snr y, then kx should have snr y.
 283     Real ka = 1.5;
 284     Real kw = 1.8;
 285     auto m2m4 = boost::math::statistics::m2m4_snr_estimator(x.begin(), x.end(), ka, kw);
 286     for(size_t i = 0; i < x.size(); ++i)
 287     {
 288         x[i] *= 4096;
 289     }
 290     auto m2m4_2 = boost::math::statistics::m2m4_snr_estimator(x.begin(), x.end(), ka, kw);
 291     BOOST_TEST(abs(m2m4 - m2m4_2) < tol);
 292 }
 293
 294 int main()
 295 {
 296     test_absolute_gini_coefficient<float>();
 297     test_absolute_gini_coefficient<double>();
 298     test_absolute_gini_coefficient<long double>();
 299
 300     test_hoyer_sparsity<float>();
 301     test_hoyer_sparsity<double>();
 302     test_hoyer_sparsity<long double>();
 303     test_hoyer_sparsity<cpp_bin_float_50>();
 304
 305     test_integer_hoyer_sparsity<int>();
 306     test_integer_hoyer_sparsity<unsigned>();
 307
 308     test_complex_hoyer_sparsity<std::complex<float>>();
 309     test_complex_hoyer_sparsity<std::complex<double>>();
 310     test_complex_hoyer_sparsity<std::complex<long double>>();
 311     test_complex_hoyer_sparsity<cpp_complex_50>();
 312
 313     test_oracle_snr<float>();
 314     test_oracle_snr<double>();
 315     test_oracle_snr<long double>();
 316     test_oracle_snr<cpp_bin_float_50>();
 317
 318     test_integer_oracle_snr<int>();
 319     test_integer_oracle_snr<unsigned>();
 320
 321     test_complex_oracle_snr<std::complex<float>>();
 322     test_complex_oracle_snr<std::complex<double>>();
 323     test_complex_oracle_snr<std::complex<long double>>();
 324     test_complex_oracle_snr<cpp_complex_50>();
 325
 326     test_m2m4_snr_estimator<float>();
 327     test_m2m4_snr_estimator<double>();
 328     test_m2m4_snr_estimator<long double>();
 329
 330     return boost::report_errors();
 331 }