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

Generated on Tue Apr 25 01:20:39 2023 (UTC) by svndigest 0.10.4pre (dev build 1607)