test/statistics.cc

Code
Comments
Other
Rev Date Author Line
116 19 Jul 04 peter 1 // $Id$
116 19 Jul 04 peter 2
675 10 Oct 06 jari 3 /*
2119 12 Dec 09 peter 4   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
831 27 Mar 07 peter 5   Copyright (C) 2005 Peter Johansson
2119 12 Dec 09 peter 6   Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
4359 23 Aug 23 peter 7   Copyright (C) 2007 Peter Johansson
4359 23 Aug 23 peter 8   Copyright (C) 2008, 2009 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 9   Copyright (C) 2010, 2011, 2012, 2013, 2014, 2016, 2018, 2020, 2021 Peter Johansson
116 19 Jul 04 peter 10
1437 25 Aug 08 peter 11   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 12
675 10 Oct 06 jari 13   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 14   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 15   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 16   License, or (at your option) any later version.
675 10 Oct 06 jari 17
675 10 Oct 06 jari 18   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 19   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 20   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 21   General Public License for more details.
675 10 Oct 06 jari 22
675 10 Oct 06 jari 23   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 24   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 25 */
675 10 Oct 06 jari 26
2881 18 Nov 12 peter 27 #include <config.h>
2881 18 Nov 12 peter 28
1248 19 Mar 08 peter 29 #include "Suite.h"
1248 19 Mar 08 peter 30
2202 21 Feb 10 peter 31 #include "yat/classifier/Target.h"
1317 21 May 08 peter 32 #include "yat/statistics/Average.h"
675 10 Oct 06 jari 33 #include "yat/statistics/utility.h"
2202 21 Feb 10 peter 34 #include "yat/statistics/tTest.h"
1835 25 Feb 09 peter 35 #include "yat/utility/DataWeight.h"
3137 28 Nov 13 peter 36 #include "yat/utility/Matrix.h"
1404 07 Aug 08 peter 37 #include "yat/utility/MatrixWeighted.h"
1120 21 Feb 08 peter 38 #include "yat/utility/Vector.h"
675 10 Oct 06 jari 39
2202 21 Feb 10 peter 40 #include <boost/concept_archetype.hpp>
3511 22 Jul 16 peter 41 #include <boost/iterator/iterator_archetypes.hpp>
2202 21 Feb 10 peter 42
1418 19 Aug 08 peter 43 #include <cmath>
116 19 Jul 04 peter 44 #include <cstdlib>
116 19 Jul 04 peter 45 #include <iostream>
1954 07 May 09 jari 46 #include <limits>
1418 19 Aug 08 peter 47 #include <map>
1418 19 Aug 08 peter 48 #include <vector>
116 19 Jul 04 peter 49
1418 19 Aug 08 peter 50 using namespace theplu::yat;
2476 14 Apr 11 peter 51 void test_benjamini_hochberg(test::Suite&);
3236 23 May 14 peter 52 void test_benjamini_hochberg_unsorted(test::Suite&);
3745 27 Jul 18 peter 53 void test_correlation(test::Suite& suite);
3135 27 Nov 13 peter 54 void test_entropy(test::Suite&);
1835 25 Feb 09 peter 55 void test_mad(test::Suite&);
3137 28 Nov 13 peter 56 void test_mutual_information(test::Suite&);
1835 25 Feb 09 peter 57
2470 12 Apr 11 peter 58 void test_median_empty(test::Suite&);
1418 19 Aug 08 peter 59 void test_percentiler(test::Suite&);
1954 07 May 09 jari 60 void test_percentiler_nan(test::Suite&);
1418 19 Aug 08 peter 61
1418 19 Aug 08 peter 62 template<typename RandomAccessIterator>
3137 28 Nov 13 peter 63 void test_percentiler(test::Suite&, RandomAccessIterator,
1418 19 Aug 08 peter 64                       RandomAccessIterator,
1418 19 Aug 08 peter 65                       double p, double correct);
1418 19 Aug 08 peter 66
1418 19 Aug 08 peter 67 template<typename RandomAccessIterator1, typename RandomAccessIterator2>
3137 28 Nov 13 peter 68 void cmp_percentiler(test::Suite&,
1418 19 Aug 08 peter 69                      RandomAccessIterator1,
3137 28 Nov 13 peter 70                      RandomAccessIterator1,
1418 19 Aug 08 peter 71                      RandomAccessIterator2,
1418 19 Aug 08 peter 72                      RandomAccessIterator2);
1418 19 Aug 08 peter 73
1248 19 Mar 08 peter 74 int main(int argc, char* argv[])
3137 28 Nov 13 peter 75 {
1248 19 Mar 08 peter 76   test::Suite suite(argc, argv);
1248 19 Mar 08 peter 77
1120 21 Feb 08 peter 78   utility::Vector gsl_vec(10);
116 19 Jul 04 peter 79   std::vector<double> data;
511 18 Feb 06 peter 80   for (unsigned int i=0; i<10; i++){
116 19 Jul 04 peter 81     data.push_back(static_cast<double>(i));
511 18 Feb 06 peter 82     gsl_vec(i)=i;
511 18 Feb 06 peter 83   }
588 21 Jun 06 markus 84
932 05 Oct 07 peter 85   double m=statistics::median(data.begin(), data.end());
932 05 Oct 07 peter 86   double m_gsl=statistics::median(gsl_vec.begin(), gsl_vec.end());
511 18 Feb 06 peter 87   if (m!=4.5 || m!=m_gsl)
1248 19 Mar 08 peter 88     suite.add(false);
2202 21 Feb 10 peter 89   if (false) {
2202 21 Feb 10 peter 90     using statistics::median;
3870 24 Feb 20 peter 91     typedef boost::iterator_archetypes::readable_lvalue_iterator_t Access;
3511 22 Jul 16 peter 92     typedef boost::random_access_traversal_tag Traversal;
3511 22 Jul 16 peter 93     boost::iterator_archetype<double, Access, Traversal> input;
3511 22 Jul 16 peter 94     double x = median(input, input);
3511 22 Jul 16 peter 95     boost::iterator_archetype<utility::DataWeight, Access, Traversal> input2;
3511 22 Jul 16 peter 96     x = median(input2, input2);
3322 06 Oct 14 peter 97     test::avoid_compiler_warning(x);
2202 21 Feb 10 peter 98   }
1500 15 Sep 08 peter 99   statistics::percentile2(data.begin(), data.end(), 100);
1039 06 Feb 08 peter 100   data.resize(1);
1039 06 Feb 08 peter 101   statistics::median(data.begin(), data.end());
1404 07 Aug 08 peter 102   // testing percentile2
1418 19 Aug 08 peter 103   test_percentiler(suite);
588 21 Jun 06 markus 104
1954 07 May 09 jari 105   // test weighted percentiler with NaNs
1954 07 May 09 jari 106   test_percentiler_nan(suite);
1954 07 May 09 jari 107
588 21 Jun 06 markus 108   double skewness_gsl=statistics::skewness(gsl_vec);
1248 19 Mar 08 peter 109   if (!suite.equal(1-skewness_gsl, 1.0) )
1248 19 Mar 08 peter 110     suite.add(false);
588 21 Jun 06 markus 111   double kurtosis_gsl=statistics::kurtosis(gsl_vec);
1671 22 Dec 08 peter 112   suite.add(suite.equal_fix(kurtosis_gsl,-1.5616363636363637113,1e-10));
1305 14 May 08 peter 113   statistics::Average func;
1305 14 May 08 peter 114   suite.add(suite.equal(func(gsl_vec.begin(), gsl_vec.end()),4.5));
1305 14 May 08 peter 115   // easiest way to get a weighted iterator
1305 14 May 08 peter 116   classifier::MatrixLookupWeighted mlw(10,20,2.0, 1.0);
1305 14 May 08 peter 117   suite.add(suite.equal(func(mlw.begin(), mlw.end()),2.0));
2202 21 Feb 10 peter 118   // do not run compiler test
2202 21 Feb 10 peter 119   if (false) {
2202 21 Feb 10 peter 120     statistics::Average average;
2202 21 Feb 10 peter 121     double x = average(boost::input_iterator_archetype<double>(),
2202 21 Feb 10 peter 122                        boost::input_iterator_archetype<double>());
3322 06 Oct 14 peter 123     test::avoid_compiler_warning(x);
3103 02 Nov 13 peter 124     using utility::DataWeight;
3103 02 Nov 13 peter 125     x = average(boost::input_iterator_archetype_no_proxy<DataWeight>(),
3103 02 Nov 13 peter 126                 boost::input_iterator_archetype_no_proxy<DataWeight>());
3322 06 Oct 14 peter 127     test::avoid_compiler_warning(x);
2202 21 Feb 10 peter 128   }
3103 02 Nov 13 peter 129
1835 25 Feb 09 peter 130   test_mad(suite);
1305 14 May 08 peter 131
2202 21 Feb 10 peter 132   // do not run compiler test
2202 21 Feb 10 peter 133   if (false) {
2202 21 Feb 10 peter 134     statistics::tTest t_test;
2202 21 Feb 10 peter 135     classifier::Target target;
3135 27 Nov 13 peter 136     add(t_test, boost::forward_iterator_archetype<double>(),
2202 21 Feb 10 peter 137         boost::forward_iterator_archetype<double>(), target);
3135 27 Nov 13 peter 138     add(t_test, boost::forward_iterator_archetype<utility::DataWeight>(),
2202 21 Feb 10 peter 139         boost::forward_iterator_archetype<utility::DataWeight>(), target);
2202 21 Feb 10 peter 140   }
2476 14 Apr 11 peter 141   test_benjamini_hochberg(suite);
3236 23 May 14 peter 142   test_benjamini_hochberg_unsorted(suite);
3135 27 Nov 13 peter 143   test_entropy(suite);
2470 12 Apr 11 peter 144   test_median_empty(suite);
3137 28 Nov 13 peter 145   test_mutual_information(suite);
3745 27 Jul 18 peter 146   test_correlation(suite);
1305 14 May 08 peter 147   return suite.return_value();
116 19 Jul 04 peter 148 }
1418 19 Aug 08 peter 149
1835 25 Feb 09 peter 150
3745 27 Jul 18 peter 151 void test_correlation(test::Suite& suite)
3745 27 Jul 18 peter 152 {
4053 26 Mar 21 peter 153   utility::Matrix x(100,10);
4053 26 Mar 21 peter 154   for (size_t i=0; i<x.rows(); ++i)
4053 26 Mar 21 peter 155     for (size_t j=0; j<x.columns(); ++j)
4053 26 Mar 21 peter 156       x(i, j) = sin(5*i + j);
3745 27 Jul 18 peter 157   utility::Matrix C = statistics::correlation(x);
3745 27 Jul 18 peter 158   if (!suite.equal(C(0,0), 1.0)) {
3745 27 Jul 18 peter 159     suite.add(false);
3745 27 Jul 18 peter 160     suite.err() << "correlation element(0, 0) not unity\n";
3745 27 Jul 18 peter 161   }
4053 26 Mar 21 peter 162
4053 26 Mar 21 peter 163   for (size_t i=0; i<C.rows(); ++i)
4053 26 Mar 21 peter 164     for (size_t j=0; j<C.columns(); ++j) {
4053 26 Mar 21 peter 165       statistics::AveragerPair ap;
4053 26 Mar 21 peter 166       add(ap, x.begin_column(i), x.end_column(i), x.begin_column(j));
4053 26 Mar 21 peter 167       if (!suite.equal(C(i, j), ap.correlation(), 10)) {
4053 26 Mar 21 peter 168         suite.add(false);
4053 26 Mar 21 peter 169         suite.err() << "error: C(" << i << ", " << j << ") is " << C(i,j)
4053 26 Mar 21 peter 170                     << "; expected " << ap.correlation() << "\n";
4053 26 Mar 21 peter 171       }
4053 26 Mar 21 peter 172     }
4053 26 Mar 21 peter 173
3745 27 Jul 18 peter 174 }
3745 27 Jul 18 peter 175
3745 27 Jul 18 peter 176
2476 14 Apr 11 peter 177 void test_benjamini_hochberg(test::Suite& suite)
2476 14 Apr 11 peter 178 {
2476 14 Apr 11 peter 179   std::vector<double> p;
2476 14 Apr 11 peter 180   p.push_back(0.0001);
2476 14 Apr 11 peter 181   p.push_back(0.01);
2476 14 Apr 11 peter 182   p.push_back(0.015);
2476 14 Apr 11 peter 183   p.push_back(0.5);
2476 14 Apr 11 peter 184   p.push_back(0.99);
2476 14 Apr 11 peter 185   std::vector<double> q(p.size());
2476 14 Apr 11 peter 186   statistics::benjamini_hochberg(p.begin(), p.end(), q.begin());
2476 14 Apr 11 peter 187   suite.add(suite.equal(q[0], p[0]*5));
2476 14 Apr 11 peter 188   suite.add(suite.equal(q[1], p[1]*2.5));
2476 14 Apr 11 peter 189   suite.add(suite.equal(q[2], 0.025));
2476 14 Apr 11 peter 190   suite.add(suite.equal(q[3], p[3]*1.25));
2476 14 Apr 11 peter 191   suite.add(suite.equal(q[4], 0.99));
2476 14 Apr 11 peter 192
2476 14 Apr 11 peter 193   // do nut run compiler test
2476 14 Apr 11 peter 194   if (false) {
2476 14 Apr 11 peter 195     using statistics::benjamini_hochberg;
3511 22 Jul 16 peter 196
3511 22 Jul 16 peter 197     typedef double Value;
3511 22 Jul 16 peter 198     typedef boost::iterator_archetypes::readable_iterator_t Access;
3511 22 Jul 16 peter 199     typedef boost::bidirectional_traversal_tag Traversal;
3511 22 Jul 16 peter 200     typedef boost::iterator_archetypes::readable_writable_iterator_t Access2;
3511 22 Jul 16 peter 201     typedef boost::bidirectional_traversal_tag Traversal2;
3511 22 Jul 16 peter 202
3511 22 Jul 16 peter 203     boost::iterator_archetype<Value, Access, Traversal> iterator;
3511 22 Jul 16 peter 204     boost::iterator_archetype<Value, Access2, Traversal2> iterator2;
3511 22 Jul 16 peter 205
3511 22 Jul 16 peter 206     benjamini_hochberg(iterator, iterator, iterator2);
2476 14 Apr 11 peter 207   }
2476 14 Apr 11 peter 208 }
2476 14 Apr 11 peter 209
2476 14 Apr 11 peter 210
3236 23 May 14 peter 211 void test_benjamini_hochberg_unsorted(test::Suite& suite)
3236 23 May 14 peter 212 {
3236 23 May 14 peter 213   std::vector<double> p;
3236 23 May 14 peter 214   p.push_back(0.015);
3236 23 May 14 peter 215   p.push_back(0.0001);
3236 23 May 14 peter 216   p.push_back(0.01);
3236 23 May 14 peter 217   p.push_back(0.5);
3236 23 May 14 peter 218   p.push_back(0.99);
3236 23 May 14 peter 219   std::vector<double> q(p.size());
3236 23 May 14 peter 220   statistics::benjamini_hochberg_unsorted(p.begin(), p.end(), q.begin());
3236 23 May 14 peter 221   suite.add(suite.equal(q[1], p[1]*5));
3236 23 May 14 peter 222   suite.add(suite.equal(q[2], p[2]*2.5));
3236 23 May 14 peter 223   suite.add(suite.equal(q[0], 0.025));
3236 23 May 14 peter 224   suite.add(suite.equal(q[3], p[3]*1.25));
3236 23 May 14 peter 225   suite.add(suite.equal(q[4], 0.99));
3236 23 May 14 peter 226
3236 23 May 14 peter 227   // do nut run compiler test
3236 23 May 14 peter 228   if (false) {
3236 23 May 14 peter 229     using statistics::benjamini_hochberg_unsorted;
3511 22 Jul 16 peter 230
3511 22 Jul 16 peter 231     typedef double Value;
3511 22 Jul 16 peter 232     typedef boost::iterator_archetypes::readable_iterator_t Access;
3511 22 Jul 16 peter 233     typedef boost::random_access_traversal_tag Traversal;
3511 22 Jul 16 peter 234     typedef boost::iterator_archetypes::readable_writable_iterator_t Access2;
3511 22 Jul 16 peter 235     typedef boost::random_access_traversal_tag Traversal2;
3511 22 Jul 16 peter 236
3511 22 Jul 16 peter 237     boost::iterator_archetype<Value, Access, Traversal> input;
3511 22 Jul 16 peter 238     boost::iterator_archetype<Value, Access2, Traversal2> result;
3511 22 Jul 16 peter 239
3236 23 May 14 peter 240     benjamini_hochberg_unsorted(input, input, result);
3236 23 May 14 peter 241   }
3236 23 May 14 peter 242 }
3236 23 May 14 peter 243
3236 23 May 14 peter 244
3135 27 Nov 13 peter 245 void test_entropy(test::Suite& suite)
3135 27 Nov 13 peter 246 {
3135 27 Nov 13 peter 247   suite.out() << "testing entropy(2)\n";
3135 27 Nov 13 peter 248   using statistics::entropy;
3135 27 Nov 13 peter 249   std::vector<int> x(10000,0);
3135 27 Nov 13 peter 250   x[512] = 42;
3135 27 Nov 13 peter 251   double e = entropy(x.begin(), x.end());
3135 27 Nov 13 peter 252   if (e>1e-15) {
3135 27 Nov 13 peter 253     suite.add(false);
3135 27 Nov 13 peter 254     suite.out() << "entropy: " << e << " expected close to 0\n";
3135 27 Nov 13 peter 255   }
3135 27 Nov 13 peter 256   x[0] = 42;
3135 27 Nov 13 peter 257   e = entropy(x.begin(), x.end());
3135 27 Nov 13 peter 258   if (e<=0) {
3135 27 Nov 13 peter 259     suite.add(false);
3135 27 Nov 13 peter 260     suite.out() << "entropy: " << e << " expected > 0\n";
3135 27 Nov 13 peter 261   }
3135 27 Nov 13 peter 262
3135 27 Nov 13 peter 263   // do not run compiler test
3135 27 Nov 13 peter 264   if (false) {
3135 27 Nov 13 peter 265     entropy(boost::input_iterator_archetype<double>(),
3135 27 Nov 13 peter 266             boost::input_iterator_archetype<double>());
3135 27 Nov 13 peter 267   }
3135 27 Nov 13 peter 268 }
3135 27 Nov 13 peter 269
3135 27 Nov 13 peter 270
1835 25 Feb 09 peter 271 void test_mad(test::Suite& suite)
1835 25 Feb 09 peter 272 {
1835 25 Feb 09 peter 273   suite.err() << "testing mad" << std::endl;
1835 25 Feb 09 peter 274   utility::Vector x(3);
1835 25 Feb 09 peter 275   x(0) = 3;
1835 25 Feb 09 peter 276   x(1) = 1;
1835 25 Feb 09 peter 277   x(2) = 100;
1835 25 Feb 09 peter 278   suite.add(suite.equal(statistics::mad(x.begin(), x.end()), 2));
1835 25 Feb 09 peter 279
1835 25 Feb 09 peter 280   std::vector<utility::DataWeight> wx(3);
1835 25 Feb 09 peter 281   wx[0] = utility::DataWeight(3, 0.4);
1835 25 Feb 09 peter 282   wx[1] = utility::DataWeight(1, 0.4);
1835 25 Feb 09 peter 283   wx[2] = utility::DataWeight(100, 0.6);
1835 25 Feb 09 peter 284   suite.add(suite.equal(statistics::mad(wx.begin(), wx.end()), 2));
2202 21 Feb 10 peter 285   // do not run compiler test
2202 21 Feb 10 peter 286   if (false) {
2202 21 Feb 10 peter 287     using statistics::mad;
3870 24 Feb 20 peter 288     typedef boost::iterator_archetypes::readable_lvalue_iterator_t Access;
3511 22 Jul 16 peter 289     typedef boost::random_access_traversal_tag Traversal;
3511 22 Jul 16 peter 290     boost::iterator_archetype<double, Access, Traversal> input;
3511 22 Jul 16 peter 291     double x = mad(input, input);
3870 24 Feb 20 peter 292     typedef boost::iterator_archetypes::readable_iterator_t Access2;
3870 24 Feb 20 peter 293     boost::iterator_archetype<utility::DataWeight, Access2, Traversal> input2;
3511 22 Jul 16 peter 294     x = mad(input2, input2);
3322 06 Oct 14 peter 295     test::avoid_compiler_warning(x);
2202 21 Feb 10 peter 296   }
1835 25 Feb 09 peter 297 }
1835 25 Feb 09 peter 298
1835 25 Feb 09 peter 299
3137 28 Nov 13 peter 300 void test_mutual_information(test::Suite& suite)
3137 28 Nov 13 peter 301 {
3137 28 Nov 13 peter 302   suite.out() << "testing mutual_information\n";
3137 28 Nov 13 peter 303   using statistics::mutual_information;
3137 28 Nov 13 peter 304   utility::Matrix x(2,2);
3137 28 Nov 13 peter 305   x(0,0) = 100;
3137 28 Nov 13 peter 306   x(1,1) = 100;
3137 28 Nov 13 peter 307   double mi = mutual_information(x);
3137 28 Nov 13 peter 308   if (!suite.add(suite.equal(mi,1.0,100))) {
3137 28 Nov 13 peter 309     suite.err() << "error: mutual information: " << mi << "\n";
3137 28 Nov 13 peter 310   }
3137 28 Nov 13 peter 311
3137 28 Nov 13 peter 312   // testing a non-square Matrix
3137 28 Nov 13 peter 313   x.resize(3,4,0);
3137 28 Nov 13 peter 314   x(0,0) = 1;
3137 28 Nov 13 peter 315   x(1,1) = 1;
3137 28 Nov 13 peter 316   x(2,2) = 1;
3137 28 Nov 13 peter 317   x(2,3) = 1;
3137 28 Nov 13 peter 318   mi = mutual_information(x);
3137 28 Nov 13 peter 319   suite.out() << "mi: " << mi << "\n";
3137 28 Nov 13 peter 320 }
3137 28 Nov 13 peter 321
3137 28 Nov 13 peter 322
2470 12 Apr 11 peter 323 // test for ticket #660
2470 12 Apr 11 peter 324 void test_median_empty(test::Suite& suite)
2470 12 Apr 11 peter 325 {
2470 12 Apr 11 peter 326   std::vector<double> x;
2470 12 Apr 11 peter 327   double m = 0;
2470 12 Apr 11 peter 328   m = statistics::median(x.begin(), x.end(), true);
3322 06 Oct 14 peter 329   test::avoid_compiler_warning(m);
2470 12 Apr 11 peter 330 }
2470 12 Apr 11 peter 331
2470 12 Apr 11 peter 332
1418 19 Aug 08 peter 333 void test_percentiler(test::Suite& suite)
1418 19 Aug 08 peter 334 {
1418 19 Aug 08 peter 335   suite.err() << "testing unweighted percentile2" << std::endl;
1418 19 Aug 08 peter 336   std::vector<double> x;
1418 19 Aug 08 peter 337   x.reserve(6);
1418 19 Aug 08 peter 338   for (unsigned int i=0; i<5; i++){
1420 20 Aug 08 peter 339     x.push_back(static_cast<double>(i+1));
1418 19 Aug 08 peter 340   }
1420 20 Aug 08 peter 341   test_percentiler(suite, x.begin(), x.end(), 50, 3);
1420 20 Aug 08 peter 342   x.push_back(6);
1420 20 Aug 08 peter 343   test_percentiler(suite, x.begin(), x.end(), 50, 3.5);
1420 20 Aug 08 peter 344   test_percentiler(suite, x.begin(), x.end(), 25, 2);
1420 20 Aug 08 peter 345   test_percentiler(suite, x.begin(), x.end(), 0, 1);
1420 20 Aug 08 peter 346   test_percentiler(suite, x.begin(), x.end(), 10, 1);
1418 19 Aug 08 peter 347
1418 19 Aug 08 peter 348   suite.err() << "testing duplication of data\n";
1418 19 Aug 08 peter 349   std::vector<double> x2(x);
1418 19 Aug 08 peter 350   for (size_t i=0; i<x.size(); ++i)
1418 19 Aug 08 peter 351     x2.push_back(x[i]);
1418 19 Aug 08 peter 352   cmp_percentiler(suite, x.begin(), x.end(), x2.begin(), x2.end());
1418 19 Aug 08 peter 353
1418 19 Aug 08 peter 354
1418 19 Aug 08 peter 355   // testing weighted
1418 19 Aug 08 peter 356
1418 19 Aug 08 peter 357   suite.err() << "testing weighted percentile2" << std::endl;
1418 19 Aug 08 peter 358   std::vector<utility::DataWeight> xw(x.size());
1418 19 Aug 08 peter 359   for (size_t i=0; i<xw.size(); ++i) {
1418 19 Aug 08 peter 360     xw[i].data() = x[i];
1418 19 Aug 08 peter 361     xw[i].weight() = 1.0;
1418 19 Aug 08 peter 362   }
1418 19 Aug 08 peter 363   const std::vector<utility::DataWeight> xw_orig(xw);
1420 20 Aug 08 peter 364   suite.err() << "testing weighted" << std::endl;
1420 20 Aug 08 peter 365   test_percentiler(suite, xw.begin(), xw.end(), 0, 1);
1420 20 Aug 08 peter 366   test_percentiler(suite, xw.begin(), xw.end(), 100, 6);
1420 20 Aug 08 peter 367   test_percentiler(suite, xw.begin(), xw.end(), 49, 3);
1420 20 Aug 08 peter 368   test_percentiler(suite, xw.begin(), xw.end(), 51, 4);
1420 20 Aug 08 peter 369   test_percentiler(suite, xw.begin(), xw.end(), 50, 3.5);
1420 20 Aug 08 peter 370   test_percentiler(suite, x.begin(), x.end(), 10, 1);
1420 20 Aug 08 peter 371
1418 19 Aug 08 peter 372   suite.err() << "testing weighted with unity weights" << std::endl;
1418 19 Aug 08 peter 373   cmp_percentiler(suite, x.begin(), x.end(), xw.begin(), xw.end());
1418 19 Aug 08 peter 374
1418 19 Aug 08 peter 375   suite.err() << "testing that w=0 equals removed data point\n";
1418 19 Aug 08 peter 376   xw=xw_orig;
1418 19 Aug 08 peter 377   std::vector<utility::DataWeight> xw2(xw_orig);
1418 19 Aug 08 peter 378   xw[3].weight() = 0.0;
1418 19 Aug 08 peter 379   xw2.erase(xw2.begin()+3);
1418 19 Aug 08 peter 380   cmp_percentiler(suite, xw.begin(), xw.end(), xw2.begin(), xw2.end());
1418 19 Aug 08 peter 381
1418 19 Aug 08 peter 382   suite.err() << "testing rescaling of weights\n";
1418 19 Aug 08 peter 383   xw2 = xw;
1418 19 Aug 08 peter 384   for (size_t i=0; i<xw2.size(); ++i)
1418 19 Aug 08 peter 385     xw2[i].weight()*=2;
1418 19 Aug 08 peter 386   cmp_percentiler(suite, xw.begin(), xw.end(), xw2.begin(), xw2.end());
1418 19 Aug 08 peter 387
2202 21 Feb 10 peter 388   // do not run compiler test
2202 21 Feb 10 peter 389   if (false) {
2202 21 Feb 10 peter 390     statistics::Percentiler percentiler(50);
3870 24 Feb 20 peter 391     typedef boost::iterator_archetypes::readable_lvalue_iterator_t Access;
3511 22 Jul 16 peter 392     typedef boost::random_access_traversal_tag Traversal;
3511 22 Jul 16 peter 393     boost::iterator_archetype<double, Access, Traversal> input;
3511 22 Jul 16 peter 394     double x = percentiler(input, input);
3511 22 Jul 16 peter 395     boost::iterator_archetype<utility::DataWeight, Access, Traversal> input2;
3511 22 Jul 16 peter 396     x = percentiler(input2, input2);
3322 06 Oct 14 peter 397     test::avoid_compiler_warning(x);
2202 21 Feb 10 peter 398   }
1418 19 Aug 08 peter 399 }
1418 19 Aug 08 peter 400
1954 07 May 09 jari 401 void test_percentiler_nan(test::Suite& suite)
1954 07 May 09 jari 402 {
1954 07 May 09 jari 403   using utility::DataWeight;
1954 07 May 09 jari 404   std::vector<double> v;
1954 07 May 09 jari 405   v.push_back(1);
1954 07 May 09 jari 406   v.push_back(10);
1954 07 May 09 jari 407   v.push_back(4);
1954 07 May 09 jari 408   v.push_back(2);
1954 07 May 09 jari 409   std::vector<DataWeight> wv(5);
1954 07 May 09 jari 410   wv[0] = DataWeight(v[0]);
1954 07 May 09 jari 411   wv[1] = DataWeight(v[1]);
1954 07 May 09 jari 412   wv[2] = DataWeight(std::numeric_limits<double>::quiet_NaN(), 0.0);
1954 07 May 09 jari 413   wv[3] = DataWeight(v[2]);
1954 07 May 09 jari 414   wv[4] = DataWeight(v[3]);
1954 07 May 09 jari 415
1954 07 May 09 jari 416   cmp_percentiler(suite, v.begin(), v.end(), wv.begin(), wv.end());
1954 07 May 09 jari 417 }
1954 07 May 09 jari 418
1418 19 Aug 08 peter 419 template<typename RandomAccessIterator>
4200 19 Aug 22 peter 420 void test_percentiler(test::Suite& suite,
4200 19 Aug 22 peter 421                       RandomAccessIterator first,
1418 19 Aug 08 peter 422                       RandomAccessIterator last,
1418 19 Aug 08 peter 423                       double p, double correct)
1418 19 Aug 08 peter 424 {
1418 19 Aug 08 peter 425   using statistics::percentile2;
1418 19 Aug 08 peter 426   double x = percentile2(first, last, p);
1418 19 Aug 08 peter 427   if (!suite.add(suite.equal(x, correct, 10))) {
1420 20 Aug 08 peter 428     suite.err() << "Error in percentile2 for " << p << "th percentile \n";
1418 19 Aug 08 peter 429     suite.err() << "  calculated value: " << x << "\n";
1418 19 Aug 08 peter 430     suite.err() << "  expected value: " << correct << "\n";
1418 19 Aug 08 peter 431   }
1418 19 Aug 08 peter 432 }
1418 19 Aug 08 peter 433
1418 19 Aug 08 peter 434 template<typename RandomAccessIterator1, typename RandomAccessIterator2>
4200 19 Aug 22 peter 435 void cmp_percentiler(test::Suite& suite,
4200 19 Aug 22 peter 436                      RandomAccessIterator1 first1,
1418 19 Aug 08 peter 437                      RandomAccessIterator1 last1,
1418 19 Aug 08 peter 438                      RandomAccessIterator2 first2,
1418 19 Aug 08 peter 439                      RandomAccessIterator2 last2)
1418 19 Aug 08 peter 440 {
2470 12 Apr 11 peter 441   for (double p=0; p<=100; p+=10) {
1418 19 Aug 08 peter 442     double correct=statistics::percentile2(first1, last1, p);
1418 19 Aug 08 peter 443     test_percentiler(suite, first2, last2, p, correct);
1418 19 Aug 08 peter 444   }
1418 19 Aug 08 peter 445
1418 19 Aug 08 peter 446 }

Generated on Sat Oct 14 02:21:37 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)