htswanalysis/MACS/lib/gsl/gsl-1.11/statistics/test_float_source.c

   1 /* statistics/test_float_source.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jim Davies, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 void FUNCTION (test, func) (const size_t stridea, const size_t strideb);
  21
  22 void
  23 FUNCTION (test, func) (const size_t stridea, const size_t strideb)
  24 {
  25   /* sample sets of doubles */
  26   size_t i;
  27   const size_t na = 14, nb = 14;
  28
  29   const double rawa[] =
  30   {.0421, .0941, .1064, .0242, .1331,
  31    .0773, .0243, .0815, .1186, .0356,
  32    .0728, .0999, .0614, .0479};
  33
  34   const double rawb[] =
  35   {.1081, .0986, .1566, .1961, .1125,
  36    .1942, .1079, .1021, .1583, .1673,
  37    .1675, .1856, .1688, .1512};
  38
  39   const double raww[] =
  40   {.0000, .0000, .0000, 3.000, .0000,
  41    1.000, 1.000, 1.000, 0.000, .5000,
  42    7.000, 5.000, 4.000, 0.123};
  43
  44   BASE * sorted ;
  45
  46   BASE * groupa = (BASE *) malloc (stridea * na * sizeof(BASE));
  47   BASE * groupb = (BASE *) malloc (strideb * nb * sizeof(BASE));
  48   BASE * w = (BASE *) malloc (strideb * na * sizeof(BASE));
  49
  50 #ifdef BASE_FLOAT
  51   double rel = 1e-6;
  52 #else
  53   double rel = 1e-10;
  54 #endif
  55
  56   for (i = 0 ; i < na ; i++)
  57     groupa[i * stridea] = (BASE) rawa[i] ;
  58
  59   for (i = 0 ; i < na ; i++)
  60     w[i * strideb] = (BASE) raww[i] ;
  61
  62   for (i = 0 ; i < nb ; i++)
  63     groupb[i * strideb] = (BASE) rawb[i] ;
  64
  65
  66   {
  67     double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
  68     double expected = 0.0728;
  69     gsl_test_rel (mean, expected, rel, NAME(gsl_stats) "_mean");
  70   }
  71
  72   {
  73     double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
  74     double var = FUNCTION(gsl_stats,variance_with_fixed_mean) (groupa, stridea, na, mean);
  75     double expected = 0.00113837428571429;
  76     gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_variance_with_fixed_mean");
  77   }
  78
  79
  80   {
  81     double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
  82     double var = FUNCTION(gsl_stats,sd_with_fixed_mean) (groupa, stridea, na, mean);
  83     double expected = 0.0337398026922845;
  84     gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_sd_with_fixed_mean");
  85   }
  86
  87
  88   {
  89     double var = FUNCTION(gsl_stats,variance) (groupb, strideb, nb);
  90     double expected = 0.00124956615384615;
  91     gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_variance");
  92   }
  93
  94   {
  95     double sd = FUNCTION(gsl_stats,sd) (groupa, stridea, na);
  96     double expected = 0.0350134479659107;
  97     gsl_test_rel (sd, expected, rel, NAME(gsl_stats) "_sd");
  98   }
  99
 100   {
 101     double ss = FUNCTION(gsl_stats,tss) (groupb, strideb, nb);
 102     double expected = 0.01624436;
 103     gsl_test_rel (ss, expected, rel, NAME(gsl_stats) "_ss");
 104   }
 105
 106   {
 107     double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
 108     double ss = FUNCTION(gsl_stats,tss_m) (groupa, stridea, na, mean);
 109     double expected = 1.59372400000000e-02;
 110     gsl_test_rel (ss, expected, rel, NAME(gsl_stats) "_ss_m");
 111   }
 112
 113   {
 114     double absdev = FUNCTION(gsl_stats,absdev) (groupa, stridea, na);
 115     double expected = 0.0287571428571429;
 116     gsl_test_rel (absdev, expected, rel, NAME(gsl_stats) "_absdev");
 117   }
 118
 119   {
 120     double skew = FUNCTION(gsl_stats,skew) (groupa, stridea, na);
 121     double expected = 0.0954642051479004;
 122     gsl_test_rel (skew, expected, rel, NAME(gsl_stats) "_skew");
 123   }
 124
 125   {
 126     double kurt = FUNCTION(gsl_stats,kurtosis) (groupa, stridea, na);
 127     double expected = -1.38583851548909 ;
 128     gsl_test_rel (kurt, expected, rel, NAME(gsl_stats) "_kurtosis");
 129   }
 130
 131   {
 132     double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
 133     double expected = 0.0678111523670601;
 134     gsl_test_rel (wmean, expected, rel, NAME(gsl_stats) "_wmean");
 135   }
 136
 137   {
 138     double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
 139     double wvar = FUNCTION(gsl_stats,wvariance_with_fixed_mean) (w, strideb, groupa, stridea, na, wmean);
 140     double expected = 0.000615793060878654;
 141     gsl_test_rel (wvar, expected, rel, NAME(gsl_stats) "_wvariance_with_fixed_mean");
 142   }
 143
 144   {
 145     double est_wvar = FUNCTION(gsl_stats,wvariance) (w, strideb, groupa, stridea, na);
 146     double expected = 0.000769562962860317;
 147     gsl_test_rel (est_wvar, expected, rel, NAME(gsl_stats) "_wvariance");
 148   }
 149
 150   {
 151     double wsd = FUNCTION(gsl_stats,wsd) (w, strideb, groupa, stridea, na);
 152     double expected = 0.0277409978706664;
 153     gsl_test_rel (wsd, expected, rel, NAME(gsl_stats) "_wsd");
 154   }
 155
 156
 157   {
 158     double wtss = FUNCTION(gsl_stats,wtss) (w, strideb, groupa, stridea, na);
 159     double expected =  1.39310864162578e-02;
 160     gsl_test_rel (wtss, expected, rel, NAME(gsl_stats) "_wtss");
 161   }
 162
 163   {
 164     double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
 165     double wtss = FUNCTION(gsl_stats,wtss_m) (w, strideb, groupa, stridea, na, wmean);
 166     double expected =  1.39310864162578e-02;
 167     gsl_test_rel (wtss, expected, rel, NAME(gsl_stats) "_wtss_m");
 168   }
 169
 170
 171
 172   {
 173     double wabsdev = FUNCTION(gsl_stats,wabsdev) (w, strideb, groupa, stridea, na);
 174     double expected = 0.0193205027504008;
 175     gsl_test_rel (wabsdev, expected, rel, NAME(gsl_stats) "_wabsdev");
 176   }
 177
 178   {
 179     double wskew = FUNCTION(gsl_stats,wskew) (w, strideb, groupa, stridea, na);
 180     double expected = -0.373631000307076;
 181     gsl_test_rel (wskew, expected, rel, NAME(gsl_stats) "_wskew");
 182   }
 183
 184   {
 185     double wkurt = FUNCTION(gsl_stats,wkurtosis) (w, strideb, groupa, stridea, na);
 186     double expected = -1.48114233353963;
 187     gsl_test_rel (wkurt, expected, rel, NAME(gsl_stats) "_wkurtosis");
 188   }
 189
 190   {
 191     double c = FUNCTION(gsl_stats,covariance) (groupa, stridea, groupb, strideb, nb);
 192     double expected = -0.000139021538461539;
 193     gsl_test_rel (c, expected, rel, NAME(gsl_stats) "_covariance");
 194   }
 195
 196   {
 197     double r = FUNCTION(gsl_stats,correlation) (groupa, stridea, groupb, strideb, nb);
 198     double expected = -0.112322712666074171;
 199     gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_correlation");
 200   }
 201
 202   {
 203     double pv = FUNCTION(gsl_stats,pvariance) (groupa, stridea, na, groupb, strideb, nb);
 204     double expected = 0.00123775384615385;
 205     gsl_test_rel (pv, expected, rel, NAME(gsl_stats) "_pvariance");
 206   }
 207
 208   {
 209     double t = FUNCTION(gsl_stats,ttest) (groupa, stridea, na, groupb, strideb, nb);
 210     double expected = -5.67026326985851;
 211     gsl_test_rel (t, expected, rel, NAME(gsl_stats) "_ttest");
 212   }
 213
 214   {
 215     BASE expected = (BASE)0.1331;
 216     gsl_test  (FUNCTION(gsl_stats,max) (groupa, stridea, na) != expected,
 217                NAME(gsl_stats) "_max (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 218                FUNCTION(gsl_stats,max) (groupa, stridea, na), expected);
 219   }
 220
 221   {
 222     BASE min = FUNCTION(gsl_stats,min) (groupa, stridea, na);
 223     BASE expected = (BASE)0.0242;
 224     gsl_test (min != expected,
 225               NAME(gsl_stats) "_min (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 226               min, expected);
 227   }
 228
 229   {
 230     BASE min, max;
 231     BASE expected_max = (BASE)0.1331;
 232     BASE expected_min = (BASE)0.0242;
 233
 234     FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na);
 235
 236     gsl_test  (max != expected_max,
 237                NAME(gsl_stats) "_minmax max (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 238                max, expected_max);
 239     gsl_test  (min != expected_min,
 240                NAME(gsl_stats) "_minmax min (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 241                min, expected_min);
 242   }
 243
 244   {
 245     int max_index = FUNCTION(gsl_stats,max_index) (groupa, stridea, na);
 246     int expected = 4;
 247     gsl_test (max_index != expected,
 248               NAME(gsl_stats) "_max_index (%d observed vs %d expected)",
 249               max_index, expected);
 250   }
 251
 252   {
 253     int min_index = FUNCTION(gsl_stats,min_index) (groupa, stridea, na);
 254     int expected = 3;
 255     gsl_test (min_index != expected,
 256               NAME(gsl_stats) "_min_index (%d observed vs %d expected)",
 257               min_index, expected);
 258   }
 259
 260   {
 261     size_t min_index, max_index;
 262     size_t expected_max_index = 4;
 263     size_t expected_min_index = 3;
 264
 265     FUNCTION(gsl_stats,minmax_index) (&min_index, &max_index, groupa, stridea, na);
 266
 267     gsl_test  (max_index != expected_max_index,
 268                NAME(gsl_stats) "_minmax_index max (%u observed vs %u expected)",
 269                max_index, expected_max_index);
 270     gsl_test  (min_index != expected_min_index,
 271                NAME(gsl_stats) "_minmax_index min (%u observed vs %u expected)",
 272                min_index, expected_min_index);
 273   }
 274
 275
 276   sorted = (BASE *) malloc(stridea * na * sizeof(BASE)) ;
 277
 278   for (i = 0 ; i < na ; i++)
 279     sorted[stridea * i] = groupa[stridea * i] ;
 280
 281   TYPE(gsl_sort)(sorted, stridea, na) ;
 282
 283   {
 284     double median = FUNCTION(gsl_stats,median_from_sorted_data)(sorted, stridea, na) ;
 285     double expected = 0.07505;
 286     gsl_test_rel  (median,expected, rel,
 287                    NAME(gsl_stats) "_median_from_sorted_data (even)");
 288   }
 289
 290   {
 291     double median = FUNCTION(gsl_stats,median_from_sorted_data)(sorted, stridea, na - 1) ;
 292     double expected = 0.0728;
 293     gsl_test_rel  (median,expected, rel,
 294                    NAME(gsl_stats) "_median_from_sorted_data");
 295   }
 296
 297
 298   {
 299     double zeroth = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 0.0) ;
 300     double expected = 0.0242;
 301     gsl_test_rel  (zeroth,expected, rel,
 302                    NAME(gsl_stats) "_quantile_from_sorted_data (0)");
 303   }
 304
 305   {
 306     double top = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 1.0) ;
 307     double expected = 0.1331;
 308     gsl_test_rel  (top,expected, rel,
 309                    NAME(gsl_stats) "_quantile_from_sorted_data (100)");
 310   }
 311
 312   {
 313     double median = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 0.5) ;
 314     double expected = 0.07505;
 315     gsl_test_rel  (median,expected, rel,
 316                    NAME(gsl_stats) "_quantile_from_sorted_data (50even)");
 317   }
 318
 319   {
 320     double median = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na - 1, 0.5);
 321     double expected = 0.0728;
 322     gsl_test_rel  (median,expected, rel,
 323                    NAME(gsl_stats) "_quantile_from_sorted_data (50odd)");
 324
 325   }
 326
 327   /* Test for IEEE handling - set third element to NaN */
 328
 329   groupa [3*stridea] = GSL_NAN;
 330
 331   {
 332     BASE max = FUNCTION(gsl_stats,max) (groupa, stridea, na);
 333     BASE expected = GSL_NAN;
 334     gsl_test  (!isnan(max),
 335                NAME(gsl_stats) "_max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 336                max, expected);
 337   }
 338
 339   {
 340     BASE min = FUNCTION(gsl_stats,min) (groupa, stridea, na);
 341     BASE expected = GSL_NAN;
 342     gsl_test (!isnan(min),
 343               NAME(gsl_stats) "_min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 344               min, expected);
 345   }
 346
 347   {
 348     BASE min, max;
 349     BASE expected_max = GSL_NAN;
 350     BASE expected_min = GSL_NAN;
 351
 352     FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na);
 353
 354     gsl_test  (!isnan(max),
 355                NAME(gsl_stats) "_minmax max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 356                max, expected_max);
 357     gsl_test  (!isnan(min),
 358                NAME(gsl_stats) "_minmax min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 359                min, expected_min);
 360   }
 361
 362 #ifdef FAST
 363   {
 364     BASE min, max;
 365     BASE expected_max = GSL_NAN;
 366     BASE expected_min = GSL_NAN;
 367
 368     FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na-1);
 369
 370     gsl_test  (!isnan(max),
 371                NAME(gsl_stats) "_minmax(-1) max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 372                max, expected_max);
 373     gsl_test  (!isnan(min),
 374                NAME(gsl_stats) "_minmax(-1) min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
 375                min, expected_min);
 376   }
 377 #endif
 378
 379
 380   {
 381     int max_index = FUNCTION(gsl_stats,max_index) (groupa, stridea, na);
 382     int expected = 3;
 383     gsl_test (max_index != expected,
 384               NAME(gsl_stats) "_max_index NaN (%d observed vs %d expected)",
 385               max_index, expected);
 386   }
 387
 388   {
 389     int min_index = FUNCTION(gsl_stats,min_index) (groupa, stridea, na);
 390     int expected = 3;
 391     gsl_test (min_index != expected,
 392               NAME(gsl_stats) "_min_index NaN (%d observed vs %d expected)",
 393               min_index, expected);
 394   }
 395
 396   {
 397     size_t min_index, max_index;
 398     size_t expected_max_index = 3;
 399     size_t expected_min_index = 3;
 400
 401     FUNCTION(gsl_stats,minmax_index) (&min_index, &max_index, groupa, stridea, na);
 402
 403     gsl_test  (max_index != expected_max_index,
 404                NAME(gsl_stats) "_minmax_index max NaN (%u observed vs %u expected)",
 405                max_index, expected_max_index);
 406     gsl_test  (min_index != expected_min_index,
 407                NAME(gsl_stats) "_minmax_index min NaN (%u observed vs %u expected)",
 408                min_index, expected_min_index);
 409   }
 410
 411   free (sorted);
 412   free (groupa);
 413   free (groupb);
 414   free (w);
 415
 416 }