ceph/src/boost/libs/math/test/univariate_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/multiprecision/cpp_bin_float.hpp>
  18 #include <boost/multiprecision/cpp_complex.hpp>
  19
  20 using boost::multiprecision::cpp_bin_float_50;
  21 using boost::multiprecision::cpp_complex_50;
  22
  23 /*
  24  * Test checklist:
  25  * 1) Does it work with multiprecision?
  26  * 2) Does it work with .cbegin()/.cend() if the data is not altered?
  27  * 3) Does it work with ublas and std::array? (Checking Eigen and Armadillo will make the CI system really unhappy.)
  28  * 4) Does it work with std::forward_list if a forward iterator is all that is required?
  29  * 5) Does it work with complex data if complex data is sensible?
  30  */
  31
  32  // To stress test, set global_seed = 0, global_size = huge.
  33  static const constexpr size_t global_seed = 0;
  34  static const constexpr size_t global_size = 128;
  35
  36 template<class T>
  37 std::vector<T> generate_random_vector(size_t size, size_t seed)
  38 {
  39     if (seed == 0)
  40     {
  41         std::random_device rd;
  42         seed = rd();
  43     }
  44     std::vector<T> v(size);
  45
  46     std::mt19937 gen(seed);
  47
  48     if constexpr (std::is_floating_point<T>::value)
  49     {
  50         std::normal_distribution<T> dis(0, 1);
  51         for (size_t i = 0; i < v.size(); ++i)
  52         {
  53          v[i] = dis(gen);
  54         }
  55         return v;
  56     }
  57     else if constexpr (std::is_integral<T>::value)
  58     {
  59         // Rescaling by larger than 2 is UB!
  60         std::uniform_int_distribution<T> dis(std::numeric_limits<T>::lowest()/2, (std::numeric_limits<T>::max)()/2);
  61         for (size_t i = 0; i < v.size(); ++i)
  62         {
  63          v[i] = dis(gen);
  64         }
  65         return v;
  66     }
  67     else if constexpr (boost::is_complex<T>::value)
  68     {
  69         std::normal_distribution<typename T::value_type> dis(0, 1);
  70         for (size_t i = 0; i < v.size(); ++i)
  71         {
  72             v[i] = {dis(gen), dis(gen)};
  73         }
  74         return v;
  75     }
  76     else if constexpr (boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
  77     {
  78         std::normal_distribution<long double> dis(0, 1);
  79         for (size_t i = 0; i < v.size(); ++i)
  80         {
  81             v[i] = {dis(gen), dis(gen)};
  82         }
  83         return v;
  84     }
  85     else if constexpr (boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point)
  86     {
  87         std::normal_distribution<long double> dis(0, 1);
  88         for (size_t i = 0; i < v.size(); ++i)
  89         {
  90             v[i] = dis(gen);
  91         }
  92         return v;
  93     }
  94     else
  95     {
  96         BOOST_ASSERT_MSG(false, "Could not identify type for random vector generation.");
  97         return v;
  98     }
  99 }
 100
 101
 102 template<class Z>
 103 void test_integer_mean()
 104 {
 105     double tol = 100*std::numeric_limits<double>::epsilon();
 106     std::vector<Z> v{1,2,3,4,5};
 107     double mu = boost::math::statistics::mean(v);
 108     BOOST_TEST(abs(mu - 3) < tol);
 109
 110     // Work with std::array?
 111     std::array<Z, 5> w{1,2,3,4,5};
 112     mu = boost::math::statistics::mean(w);
 113     BOOST_TEST(abs(mu - 3) < tol);
 114
 115     v = generate_random_vector<Z>(global_size, global_seed);
 116     Z scale = 2;
 117
 118     double m1 = scale*boost::math::statistics::mean(v);
 119     for (auto & x : v)
 120     {
 121         x *= scale;
 122     }
 123     double m2 = boost::math::statistics::mean(v);
 124     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 125 }
 126
 127 template<class RandomAccessContainer>
 128 auto naive_mean(RandomAccessContainer const & v)
 129 {
 130     typename RandomAccessContainer::value_type sum = 0;
 131     for (auto & x : v)
 132     {
 133         sum += x;
 134     }
 135     return sum/v.size();
 136 }
 137
 138 template<class Real>
 139 void test_mean()
 140 {
 141     Real tol = std::numeric_limits<Real>::epsilon();
 142     std::vector<Real> v{1,2,3,4,5};
 143     Real mu = boost::math::statistics::mean(v.begin(), v.end());
 144     BOOST_TEST(abs(mu - 3) < tol);
 145
 146     // Does range call work?
 147     mu = boost::math::statistics::mean(v);
 148     BOOST_TEST(abs(mu - 3) < tol);
 149
 150     // Can we successfully average only part of the vector?
 151     mu = boost::math::statistics::mean(v.begin(), v.begin() + 3);
 152     BOOST_TEST(abs(mu - 2) < tol);
 153
 154     // Does it work when we const qualify?
 155     mu = boost::math::statistics::mean(v.cbegin(), v.cend());
 156     BOOST_TEST(abs(mu - 3) < tol);
 157
 158     // Does it work for std::array?
 159     std::array<Real, 7> u{1,2,3,4,5,6,7};
 160     mu = boost::math::statistics::mean(u.begin(), u.end());
 161     BOOST_TEST(abs(mu - 4) < 10*tol);
 162
 163     // Does it work for a forward iterator?
 164     std::forward_list<Real> l{1,2,3,4,5,6,7};
 165     mu = boost::math::statistics::mean(l.begin(), l.end());
 166     BOOST_TEST(abs(mu - 4) < tol);
 167
 168     // Does it work with ublas vectors?
 169     boost::numeric::ublas::vector<Real> w(7);
 170     for (size_t i = 0; i < w.size(); ++i)
 171     {
 172         w[i] = i+1;
 173     }
 174     mu = boost::math::statistics::mean(w.cbegin(), w.cend());
 175     BOOST_TEST(abs(mu - 4) < tol);
 176
 177     v = generate_random_vector<Real>(global_size, global_seed);
 178     Real scale = 2;
 179     Real m1 = scale*boost::math::statistics::mean(v);
 180     for (auto & x : v)
 181     {
 182         x *= scale;
 183     }
 184     Real m2 = boost::math::statistics::mean(v);
 185     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 186
 187     // Stress test:
 188     for (size_t i = 1; i < 30; ++i)
 189     {
 190         v = generate_random_vector<Real>(i, 12803);
 191         auto naive_ = naive_mean(v);
 192         auto higham_ = boost::math::statistics::mean(v);
 193         if (abs(higham_ - naive_) >= 100*tol*abs(naive_))
 194         {
 195             std::cout << std::hexfloat;
 196             std::cout << "Terms = " << v.size() << "\n";
 197             std::cout << "higham = " << higham_ << "\n";
 198             std::cout << "naive_ = " << naive_ << "\n";
 199         }
 200         BOOST_TEST(abs(higham_ - naive_) < 100*tol*abs(naive_));
 201     }
 202
 203 }
 204
 205 template<class Complex>
 206 void test_complex_mean()
 207 {
 208     typedef typename Complex::value_type Real;
 209     Real tol = std::numeric_limits<Real>::epsilon();
 210     std::vector<Complex> v{{0,1},{0,2},{0,3},{0,4},{0,5}};
 211     auto mu = boost::math::statistics::mean(v.begin(), v.end());
 212     BOOST_TEST(abs(mu.imag() - 3) < tol);
 213     BOOST_TEST(abs(mu.real()) < tol);
 214
 215     // Does range work?
 216     mu = boost::math::statistics::mean(v);
 217     BOOST_TEST(abs(mu.imag() - 3) < tol);
 218     BOOST_TEST(abs(mu.real()) < tol);
 219 }
 220
 221 template<class Real>
 222 void test_variance()
 223 {
 224     Real tol = std::numeric_limits<Real>::epsilon();
 225     std::vector<Real> v{1,1,1,1,1,1};
 226     Real sigma_sq = boost::math::statistics::variance(v.begin(), v.end());
 227     BOOST_TEST(abs(sigma_sq) < tol);
 228
 229     sigma_sq = boost::math::statistics::variance(v);
 230     BOOST_TEST(abs(sigma_sq) < tol);
 231
 232     Real s_sq = boost::math::statistics::sample_variance(v);
 233     BOOST_TEST(abs(s_sq) < tol);
 234
 235     std::vector<Real> u{1};
 236     sigma_sq = boost::math::statistics::variance(u.cbegin(), u.cend());
 237     BOOST_TEST(abs(sigma_sq) < tol);
 238
 239     std::array<Real, 8> w{0,1,0,1,0,1,0,1};
 240     sigma_sq = boost::math::statistics::variance(w.begin(), w.end());
 241     BOOST_TEST(abs(sigma_sq - 1.0/4.0) < tol);
 242
 243     sigma_sq = boost::math::statistics::variance(w);
 244     BOOST_TEST(abs(sigma_sq - 1.0/4.0) < tol);
 245
 246     std::forward_list<Real> l{0,1,0,1,0,1,0,1};
 247     sigma_sq = boost::math::statistics::variance(l.begin(), l.end());
 248     BOOST_TEST(abs(sigma_sq - 1.0/4.0) < tol);
 249
 250     v = generate_random_vector<Real>(global_size, global_seed);
 251     Real scale = 2;
 252     Real m1 = scale*scale*boost::math::statistics::variance(v);
 253     for (auto & x : v)
 254     {
 255         x *= scale;
 256     }
 257     Real m2 = boost::math::statistics::variance(v);
 258     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 259
 260     // Wikipedia example for a variance of N sided die:
 261     // https://en.wikipedia.org/wiki/Variance
 262     for (size_t j = 16; j < 2048; j *= 2)
 263     {
 264         v.resize(1024);
 265         Real n = v.size();
 266         for (size_t i = 0; i < v.size(); ++i)
 267         {
 268             v[i] = i + 1;
 269         }
 270
 271         sigma_sq = boost::math::statistics::variance(v);
 272
 273         BOOST_TEST(abs(sigma_sq - (n*n-1)/Real(12)) <= tol*sigma_sq);
 274     }
 275
 276 }
 277
 278 template<class Z>
 279 void test_integer_variance()
 280 {
 281     double tol = std::numeric_limits<double>::epsilon();
 282     std::vector<Z> v{1,1,1,1,1,1};
 283     double sigma_sq = boost::math::statistics::variance(v);
 284     BOOST_TEST(abs(sigma_sq) < tol);
 285
 286     std::forward_list<Z> l{0,1,0,1,0,1,0,1};
 287     sigma_sq = boost::math::statistics::variance(l.begin(), l.end());
 288     BOOST_TEST(abs(sigma_sq - 1.0/4.0) < tol);
 289
 290     v = generate_random_vector<Z>(global_size, global_seed);
 291     Z scale = 2;
 292     double m1 = scale*scale*boost::math::statistics::variance(v);
 293     for (auto & x : v)
 294     {
 295         x *= scale;
 296     }
 297     double m2 = boost::math::statistics::variance(v);
 298     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 299 }
 300
 301 template<class Z>
 302 void test_integer_skewness()
 303 {
 304     double tol = std::numeric_limits<double>::epsilon();
 305     std::vector<Z> v{1,1,1};
 306     double skew = boost::math::statistics::skewness(v);
 307     BOOST_TEST(abs(skew) < tol);
 308
 309     // Dataset is symmetric about the mean:
 310     v = {1,2,3,4,5};
 311     skew = boost::math::statistics::skewness(v);
 312     BOOST_TEST(abs(skew) < tol);
 313
 314     v = {0,0,0,0,5};
 315     // mu = 1, sigma^2 = 4, sigma = 2, skew = 3/2
 316     skew = boost::math::statistics::skewness(v);
 317     BOOST_TEST(abs(skew - 3.0/2.0) < tol);
 318
 319     std::forward_list<Z> v2{0,0,0,0,5};
 320     skew = boost::math::statistics::skewness(v);
 321     BOOST_TEST(abs(skew - 3.0/2.0) < tol);
 322
 323
 324     v = generate_random_vector<Z>(global_size, global_seed);
 325     Z scale = 2;
 326     double m1 = boost::math::statistics::skewness(v);
 327     for (auto & x : v)
 328     {
 329         x *= scale;
 330     }
 331     double m2 = boost::math::statistics::skewness(v);
 332     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 333
 334 }
 335
 336 template<class Real>
 337 void test_skewness()
 338 {
 339     Real tol = std::numeric_limits<Real>::epsilon();
 340     std::vector<Real> v{1,1,1};
 341     Real skew = boost::math::statistics::skewness(v);
 342     BOOST_TEST(abs(skew) < tol);
 343
 344     // Dataset is symmetric about the mean:
 345     v = {1,2,3,4,5};
 346     skew = boost::math::statistics::skewness(v);
 347     BOOST_TEST(abs(skew) < tol);
 348
 349     v = {0,0,0,0,5};
 350     // mu = 1, sigma^2 = 4, sigma = 2, skew = 3/2
 351     skew = boost::math::statistics::skewness(v);
 352     BOOST_TEST(abs(skew - Real(3)/Real(2)) < tol);
 353
 354     std::array<Real, 5> w1{0,0,0,0,5};
 355     skew = boost::math::statistics::skewness(w1);
 356     BOOST_TEST(abs(skew - Real(3)/Real(2)) < tol);
 357
 358     std::forward_list<Real> w2{0,0,0,0,5};
 359     skew = boost::math::statistics::skewness(w2);
 360     BOOST_TEST(abs(skew - Real(3)/Real(2)) < tol);
 361
 362     v = generate_random_vector<Real>(global_size, global_seed);
 363     Real scale = 2;
 364     Real m1 = boost::math::statistics::skewness(v);
 365     for (auto & x : v)
 366     {
 367         x *= scale;
 368     }
 369     Real m2 = boost::math::statistics::skewness(v);
 370     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 371 }
 372
 373 template<class Real>
 374 void test_kurtosis()
 375 {
 376     Real tol = std::numeric_limits<Real>::epsilon();
 377     std::vector<Real> v{1,1,1};
 378     Real kurt = boost::math::statistics::kurtosis(v);
 379     BOOST_TEST(abs(kurt) < tol);
 380
 381     v = {1,2,3,4,5};
 382     // mu =1, sigma^2 = 2, kurtosis = 17/10
 383     kurt = boost::math::statistics::kurtosis(v);
 384     BOOST_TEST(abs(kurt - Real(17)/Real(10)) < 10*tol);
 385
 386     v = {0,0,0,0,5};
 387     // mu = 1, sigma^2 = 4, sigma = 2, skew = 3/2, kurtosis = 13/4
 388     kurt = boost::math::statistics::kurtosis(v);
 389     BOOST_TEST(abs(kurt - Real(13)/Real(4)) < tol);
 390
 391     std::array<Real, 5> v1{0,0,0,0,5};
 392     kurt = boost::math::statistics::kurtosis(v1);
 393     BOOST_TEST(abs(kurt - Real(13)/Real(4)) < tol);
 394
 395     std::forward_list<Real> v2{0,0,0,0,5};
 396     kurt = boost::math::statistics::kurtosis(v2);
 397     BOOST_TEST(abs(kurt - Real(13)/Real(4)) < tol);
 398
 399     std::vector<Real> v3(10000);
 400     std::mt19937 gen(42);
 401     std::normal_distribution<long double> dis(0, 1);
 402     for (size_t i = 0; i < v3.size(); ++i) {
 403         v3[i] = dis(gen);
 404     }
 405     kurt = boost::math::statistics::kurtosis(v3);
 406     BOOST_TEST(abs(kurt - 3) < 0.1);
 407
 408     std::uniform_real_distribution<long double> udis(-1, 3);
 409     for (size_t i = 0; i < v3.size(); ++i) {
 410         v3[i] = udis(gen);
 411     }
 412     auto excess_kurtosis = boost::math::statistics::excess_kurtosis(v3);
 413     BOOST_TEST(abs(excess_kurtosis + 6.0/5.0) < 0.2);
 414
 415     v = generate_random_vector<Real>(global_size, global_seed);
 416     Real scale = 2;
 417     Real m1 = boost::math::statistics::kurtosis(v);
 418     for (auto & x : v)
 419     {
 420         x *= scale;
 421     }
 422     Real m2 = boost::math::statistics::kurtosis(v);
 423     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 424
 425     // This test only passes when there are a large number of samples.
 426     // Otherwise, the distribution doesn't generate enough outliers to give,
 427     // or generates too many, giving pretty wildly different values of kurtosis on different runs.
 428     // However, by kicking up the samples to 1,000,000, I got very close to 6 for the excess kurtosis on every run.
 429     // The CI system, however, would die on a million long doubles.
 430     //v3.resize(1000000);
 431     //std::exponential_distribution<long double> edis(0.1);
 432     //for (size_t i = 0; i < v3.size(); ++i) {
 433     //    v3[i] = edis(gen);
 434     //}
 435     //excess_kurtosis = boost::math::statistics::kurtosis(v3) - 3;
 436     //BOOST_TEST(abs(excess_kurtosis - 6.0) < 0.2);
 437 }
 438
 439 template<class Z>
 440 void test_integer_kurtosis()
 441 {
 442     double tol = std::numeric_limits<double>::epsilon();
 443     std::vector<Z> v{1,1,1};
 444     double kurt = boost::math::statistics::kurtosis(v);
 445     BOOST_TEST(abs(kurt) < tol);
 446
 447     v = {1,2,3,4,5};
 448     // mu =1, sigma^2 = 2, kurtosis = 17/10
 449     kurt = boost::math::statistics::kurtosis(v);
 450     BOOST_TEST(abs(kurt - 17.0/10.0) < 10*tol);
 451
 452     v = {0,0,0,0,5};
 453     // mu = 1, sigma^2 = 4, sigma = 2, skew = 3/2, kurtosis = 13/4
 454     kurt = boost::math::statistics::kurtosis(v);
 455     BOOST_TEST(abs(kurt - 13.0/4.0) < tol);
 456
 457     v = generate_random_vector<Z>(global_size, global_seed);
 458     Z scale = 2;
 459     double m1 = boost::math::statistics::kurtosis(v);
 460     for (auto & x : v)
 461     {
 462         x *= scale;
 463     }
 464     double m2 = boost::math::statistics::kurtosis(v);
 465     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
 466 }
 467
 468 template<class Real>
 469 void test_first_four_moments()
 470 {
 471     Real tol = 10*std::numeric_limits<Real>::epsilon();
 472     std::vector<Real> v{1,1,1};
 473     auto [M1_1, M2_1, M3_1, M4_1] = boost::math::statistics::first_four_moments(v);
 474     BOOST_TEST(abs(M1_1 - 1) < tol);
 475     BOOST_TEST(abs(M2_1) < tol);
 476     BOOST_TEST(abs(M3_1) < tol);
 477     BOOST_TEST(abs(M4_1) < tol);
 478
 479     v = {1, 2, 3, 4, 5};
 480     auto [M1_2, M2_2, M3_2, M4_2] = boost::math::statistics::first_four_moments(v);
 481     BOOST_TEST(abs(M1_2 - 3) < tol);
 482     BOOST_TEST(abs(M2_2 - 2) < tol);
 483     BOOST_TEST(abs(M3_2) < tol);
 484     BOOST_TEST(abs(M4_2 - Real(34)/Real(5)) < tol);
 485 }
 486
 487 template<class Real>
 488 void test_median()
 489 {
 490     std::mt19937 g(12);
 491     std::vector<Real> v{1,2,3,4,5,6,7};
 492
 493     Real m = boost::math::statistics::median(v.begin(), v.end());
 494     BOOST_TEST_EQ(m, 4);
 495
 496     std::shuffle(v.begin(), v.end(), g);
 497     // Does range call work?
 498     m = boost::math::statistics::median(v);
 499     BOOST_TEST_EQ(m, 4);
 500
 501     v = {1,2,3,3,4,5};
 502     m = boost::math::statistics::median(v.begin(), v.end());
 503     BOOST_TEST_EQ(m, 3);
 504     std::shuffle(v.begin(), v.end(), g);
 505     m = boost::math::statistics::median(v.begin(), v.end());
 506     BOOST_TEST_EQ(m, 3);
 507
 508     v = {1};
 509     m = boost::math::statistics::median(v.begin(), v.end());
 510     BOOST_TEST_EQ(m, 1);
 511
 512     v = {1,1};
 513     m = boost::math::statistics::median(v.begin(), v.end());
 514     BOOST_TEST_EQ(m, 1);
 515
 516     v = {2,4};
 517     m = boost::math::statistics::median(v.begin(), v.end());
 518     BOOST_TEST_EQ(m, 3);
 519
 520     v = {1,1,1};
 521     m = boost::math::statistics::median(v.begin(), v.end());
 522     BOOST_TEST_EQ(m, 1);
 523
 524     v = {1,2,3};
 525     m = boost::math::statistics::median(v.begin(), v.end());
 526     BOOST_TEST_EQ(m, 2);
 527     std::shuffle(v.begin(), v.end(), g);
 528     m = boost::math::statistics::median(v.begin(), v.end());
 529     BOOST_TEST_EQ(m, 2);
 530
 531     // Does it work with std::array?
 532     std::array<Real, 3> w{1,2,3};
 533     m = boost::math::statistics::median(w);
 534     BOOST_TEST_EQ(m, 2);
 535
 536     // Does it work with ublas?
 537     boost::numeric::ublas::vector<Real> w1(3);
 538     w1[0] = 1;
 539     w1[1] = 2;
 540     w1[2] = 3;
 541     m = boost::math::statistics::median(w);
 542     BOOST_TEST_EQ(m, 2);
 543 }
 544
 545 template<class Real>
 546 void test_median_absolute_deviation()
 547 {
 548     std::vector<Real> v{-1, 2, -3, 4, -5, 6, -7};
 549
 550     Real m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 551     BOOST_TEST_EQ(m, 4);
 552
 553     std::mt19937 g(12);
 554     std::shuffle(v.begin(), v.end(), g);
 555     m = boost::math::statistics::median_absolute_deviation(v, 0);
 556     BOOST_TEST_EQ(m, 4);
 557
 558     v = {1, -2, -3, 3, -4, -5};
 559     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 560     BOOST_TEST_EQ(m, 3);
 561     std::shuffle(v.begin(), v.end(), g);
 562     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 563     BOOST_TEST_EQ(m, 3);
 564
 565     v = {-1};
 566     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 567     BOOST_TEST_EQ(m, 1);
 568
 569     v = {-1, 1};
 570     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 571     BOOST_TEST_EQ(m, 1);
 572     // The median is zero, so coincides with the default:
 573     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end());
 574     BOOST_TEST_EQ(m, 1);
 575
 576     m = boost::math::statistics::median_absolute_deviation(v);
 577     BOOST_TEST_EQ(m, 1);
 578
 579
 580     v = {2, -4};
 581     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 582     BOOST_TEST_EQ(m, 3);
 583
 584     v = {1, -1, 1};
 585     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 586     BOOST_TEST_EQ(m, 1);
 587
 588     v = {1, 2, -3};
 589     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 590     BOOST_TEST_EQ(m, 2);
 591     std::shuffle(v.begin(), v.end(), g);
 592     m = boost::math::statistics::median_absolute_deviation(v.begin(), v.end(), 0);
 593     BOOST_TEST_EQ(m, 2);
 594
 595     std::array<Real, 3> w{1, 2, -3};
 596     m = boost::math::statistics::median_absolute_deviation(w, 0);
 597     BOOST_TEST_EQ(m, 2);
 598
 599     // boost.ublas vector?
 600     boost::numeric::ublas::vector<Real> u(6);
 601     u[0] = 1;
 602     u[1] = 2;
 603     u[2] = -3;
 604     u[3] = 1;
 605     u[4] = 2;
 606     u[5] = -3;
 607     m = boost::math::statistics::median_absolute_deviation(u, 0);
 608     BOOST_TEST_EQ(m, 2);
 609 }
 610
 611
 612 template<class Real>
 613 void test_sample_gini_coefficient()
 614 {
 615     Real tol = std::numeric_limits<Real>::epsilon();
 616     std::vector<Real> v{1,0,0};
 617     Real gini = boost::math::statistics::sample_gini_coefficient(v.begin(), v.end());
 618     BOOST_TEST(abs(gini - 1) < tol);
 619
 620     gini = boost::math::statistics::sample_gini_coefficient(v);
 621     BOOST_TEST(abs(gini - 1) < tol);
 622
 623     v[0] = 1;
 624     v[1] = 1;
 625     v[2] = 1;
 626     gini = boost::math::statistics::sample_gini_coefficient(v.begin(), v.end());
 627     BOOST_TEST(abs(gini) < tol);
 628
 629     v[0] = 0;
 630     v[1] = 0;
 631     v[2] = 0;
 632     gini = boost::math::statistics::sample_gini_coefficient(v.begin(), v.end());
 633     BOOST_TEST(abs(gini) < tol);
 634
 635     std::array<Real, 3> w{0,0,0};
 636     gini = boost::math::statistics::sample_gini_coefficient(w);
 637     BOOST_TEST(abs(gini) < tol);
 638 }
 639
 640
 641 template<class Real>
 642 void test_gini_coefficient()
 643 {
 644     Real tol = std::numeric_limits<Real>::epsilon();
 645     std::vector<Real> v{1,0,0};
 646     Real gini = boost::math::statistics::gini_coefficient(v.begin(), v.end());
 647     Real expected = Real(2)/Real(3);
 648     BOOST_TEST(abs(gini - expected) < tol);
 649
 650     gini = boost::math::statistics::gini_coefficient(v);
 651     BOOST_TEST(abs(gini - expected) < tol);
 652
 653     v[0] = 1;
 654     v[1] = 1;
 655     v[2] = 1;
 656     gini = boost::math::statistics::gini_coefficient(v.begin(), v.end());
 657     BOOST_TEST(abs(gini) < tol);
 658
 659     v[0] = 0;
 660     v[1] = 0;
 661     v[2] = 0;
 662     gini = boost::math::statistics::gini_coefficient(v.begin(), v.end());
 663     BOOST_TEST(abs(gini) < tol);
 664
 665     std::array<Real, 3> w{0,0,0};
 666     gini = boost::math::statistics::gini_coefficient(w);
 667     BOOST_TEST(abs(gini) < tol);
 668
 669     boost::numeric::ublas::vector<Real> w1(3);
 670     w1[0] = 1;
 671     w1[1] = 1;
 672     w1[2] = 1;
 673     gini = boost::math::statistics::gini_coefficient(w1);
 674     BOOST_TEST(abs(gini) < tol);
 675
 676     std::mt19937 gen(18);
 677     // Gini coefficient for a uniform distribution is (b-a)/(3*(b+a));
 678     std::uniform_real_distribution<long double> dis(0, 3);
 679     expected = (dis.b() - dis.a())/(3*(dis.b()+ dis.a()));
 680     v.resize(1024);
 681     for(size_t i = 0; i < v.size(); ++i)
 682     {
 683         v[i] = dis(gen);
 684     }
 685     gini = boost::math::statistics::gini_coefficient(v);
 686     BOOST_TEST(abs(gini - expected) < 0.01);
 687
 688 }
 689
 690 template<class Z>
 691 void test_integer_gini_coefficient()
 692 {
 693     double tol = std::numeric_limits<double>::epsilon();
 694     std::vector<Z> v{1,0,0};
 695     double gini = boost::math::statistics::gini_coefficient(v.begin(), v.end());
 696     double expected = 2.0/3.0;
 697     BOOST_TEST(abs(gini - expected) < tol);
 698
 699     gini = boost::math::statistics::gini_coefficient(v);
 700     BOOST_TEST(abs(gini - expected) < tol);
 701
 702     v[0] = 1;
 703     v[1] = 1;
 704     v[2] = 1;
 705     gini = boost::math::statistics::gini_coefficient(v.begin(), v.end());
 706     BOOST_TEST(abs(gini) < tol);
 707
 708     v[0] = 0;
 709     v[1] = 0;
 710     v[2] = 0;
 711     gini = boost::math::statistics::gini_coefficient(v.begin(), v.end());
 712     BOOST_TEST(abs(gini) < tol);
 713
 714     std::array<Z, 3> w{0,0,0};
 715     gini = boost::math::statistics::gini_coefficient(w);
 716     BOOST_TEST(abs(gini) < tol);
 717
 718     boost::numeric::ublas::vector<Z> w1(3);
 719     w1[0] = 1;
 720     w1[1] = 1;
 721     w1[2] = 1;
 722     gini = boost::math::statistics::gini_coefficient(w1);
 723     BOOST_TEST(abs(gini) < tol);
 724 }
 725
 726 template<typename Real>
 727 void test_interquartile_range()
 728 {
 729     std::mt19937 gen(486);
 730     Real iqr;
 731     // Taken from Wikipedia's example:
 732     std::vector<Real> v{7, 7, 31, 31, 47, 75, 87, 115, 116, 119, 119, 155, 177};
 733
 734     // Q1 = 31, Q3 = 119, Q3 - Q1 = 88.
 735     iqr = boost::math::statistics::interquartile_range(v);
 736     BOOST_TEST_EQ(iqr, 88);
 737
 738     std::shuffle(v.begin(), v.end(), gen);
 739     iqr = boost::math::statistics::interquartile_range(v);
 740     BOOST_TEST_EQ(iqr, 88);
 741
 742     std::shuffle(v.begin(), v.end(), gen);
 743     iqr = boost::math::statistics::interquartile_range(v);
 744     BOOST_TEST_EQ(iqr, 88);
 745
 746     std::fill(v.begin(), v.end(), 1);
 747     iqr = boost::math::statistics::interquartile_range(v);
 748     BOOST_TEST_EQ(iqr, 0);
 749
 750     v = {1,2,3};
 751     iqr = boost::math::statistics::interquartile_range(v);
 752     BOOST_TEST_EQ(iqr, 2);
 753     std::shuffle(v.begin(), v.end(), gen);
 754     iqr = boost::math::statistics::interquartile_range(v);
 755     BOOST_TEST_EQ(iqr, 2);
 756
 757     v = {0, 3, 5};
 758     iqr = boost::math::statistics::interquartile_range(v);
 759     BOOST_TEST_EQ(iqr, 5);
 760     std::shuffle(v.begin(), v.end(), gen);
 761     iqr = boost::math::statistics::interquartile_range(v);
 762     BOOST_TEST_EQ(iqr, 5);
 763
 764     v = {1,2,3,4};
 765     iqr = boost::math::statistics::interquartile_range(v);
 766     BOOST_TEST_EQ(iqr, 2);
 767     std::shuffle(v.begin(), v.end(), gen);
 768     iqr = boost::math::statistics::interquartile_range(v);
 769     BOOST_TEST_EQ(iqr, 2);
 770
 771     v = {1,2,3,4,5};
 772     // Q1 = 1.5, Q3 = 4.5
 773     iqr = boost::math::statistics::interquartile_range(v);
 774     BOOST_TEST_EQ(iqr, 3);
 775     std::shuffle(v.begin(), v.end(), gen);
 776     iqr = boost::math::statistics::interquartile_range(v);
 777     BOOST_TEST_EQ(iqr, 3);
 778
 779     v = {1,2,3,4,5,6};
 780     // Q1 = 2, Q3 = 5
 781     iqr = boost::math::statistics::interquartile_range(v);
 782     BOOST_TEST_EQ(iqr, 3);
 783     std::shuffle(v.begin(), v.end(), gen);
 784     iqr = boost::math::statistics::interquartile_range(v);
 785     BOOST_TEST_EQ(iqr, 3);
 786
 787     v = {1,2,3, 4, 5,6,7};
 788     // Q1 = 2, Q3 = 6
 789     iqr = boost::math::statistics::interquartile_range(v);
 790     BOOST_TEST_EQ(iqr, 4);
 791     std::shuffle(v.begin(), v.end(), gen);
 792     iqr = boost::math::statistics::interquartile_range(v);
 793     BOOST_TEST_EQ(iqr, 4);
 794
 795     v = {1,2,3,4,5,6,7,8};
 796     // Q1 = 2.5, Q3 = 6.5
 797     iqr = boost::math::statistics::interquartile_range(v);
 798     BOOST_TEST_EQ(iqr, 4);
 799     std::shuffle(v.begin(), v.end(), gen);
 800     iqr = boost::math::statistics::interquartile_range(v);
 801     BOOST_TEST_EQ(iqr, 4);
 802
 803     v = {1,2,3,4,5,6,7,8,9};
 804     // Q1 = 2.5, Q3 = 7.5
 805     iqr = boost::math::statistics::interquartile_range(v);
 806     BOOST_TEST_EQ(iqr, 5);
 807     std::shuffle(v.begin(), v.end(), gen);
 808     iqr = boost::math::statistics::interquartile_range(v);
 809     BOOST_TEST_EQ(iqr, 5);
 810
 811     v = {1,2,3,4,5,6,7,8,9,10};
 812     // Q1 = 3, Q3 = 8
 813     iqr = boost::math::statistics::interquartile_range(v);
 814     BOOST_TEST_EQ(iqr, 5);
 815     std::shuffle(v.begin(), v.end(), gen);
 816     iqr = boost::math::statistics::interquartile_range(v);
 817     BOOST_TEST_EQ(iqr, 5);
 818
 819     v = {1,2,3,4,5,6,7,8,9,10,11};
 820     // Q1 = 3, Q3 = 9
 821     iqr = boost::math::statistics::interquartile_range(v);
 822     BOOST_TEST_EQ(iqr, 6);
 823     std::shuffle(v.begin(), v.end(), gen);
 824     iqr = boost::math::statistics::interquartile_range(v);
 825     BOOST_TEST_EQ(iqr, 6);
 826
 827     v = {1,2,3,4,5,6,7,8,9,10,11,12};
 828     // Q1 = 3.5, Q3 = 9.5
 829     iqr = boost::math::statistics::interquartile_range(v);
 830     BOOST_TEST_EQ(iqr, 6);
 831     std::shuffle(v.begin(), v.end(), gen);
 832     iqr = boost::math::statistics::interquartile_range(v);
 833     BOOST_TEST_EQ(iqr, 6);
 834 }
 835
 836 int main()
 837 {
 838     test_mean<float>();
 839     test_mean<double>();
 840     test_mean<long double>();
 841     test_mean<cpp_bin_float_50>();
 842
 843     test_integer_mean<unsigned>();
 844     test_integer_mean<int>();
 845
 846     test_complex_mean<std::complex<float>>();
 847     test_complex_mean<cpp_complex_50>();
 848
 849     test_variance<float>();
 850     test_variance<double>();
 851     test_variance<long double>();
 852     test_variance<cpp_bin_float_50>();
 853
 854     test_integer_variance<int>();
 855     test_integer_variance<unsigned>();
 856
 857     test_skewness<float>();
 858     test_skewness<double>();
 859     test_skewness<long double>();
 860     test_skewness<cpp_bin_float_50>();
 861
 862     test_integer_skewness<int>();
 863     test_integer_skewness<unsigned>();
 864
 865     test_first_four_moments<float>();
 866     test_first_four_moments<double>();
 867     test_first_four_moments<long double>();
 868     test_first_four_moments<cpp_bin_float_50>();
 869
 870     test_kurtosis<float>();
 871     test_kurtosis<double>();
 872     test_kurtosis<long double>();
 873     // Kinda expensive:
 874     //test_kurtosis<cpp_bin_float_50>();
 875
 876     test_integer_kurtosis<int>();
 877     test_integer_kurtosis<unsigned>();
 878
 879     test_median<float>();
 880     test_median<double>();
 881     test_median<long double>();
 882     test_median<cpp_bin_float_50>();
 883     test_median<int>();
 884
 885     test_median_absolute_deviation<float>();
 886     test_median_absolute_deviation<double>();
 887     test_median_absolute_deviation<long double>();
 888     test_median_absolute_deviation<cpp_bin_float_50>();
 889
 890     test_gini_coefficient<float>();
 891     test_gini_coefficient<double>();
 892     test_gini_coefficient<long double>();
 893     test_gini_coefficient<cpp_bin_float_50>();
 894
 895     test_integer_gini_coefficient<unsigned>();
 896     test_integer_gini_coefficient<int>();
 897
 898     test_sample_gini_coefficient<float>();
 899     test_sample_gini_coefficient<double>();
 900     test_sample_gini_coefficient<long double>();
 901     test_sample_gini_coefficient<cpp_bin_float_50>();
 902
 903     test_interquartile_range<double>();
 904     test_interquartile_range<cpp_bin_float_50>();
 905     return boost::report_errors();
 906 }