/* statistics/test_float_source.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jim Davies, Brian Gough
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
void FUNCTION (test, func) (const size_t stridea, const size_t strideb);
void
FUNCTION (test, func) (const size_t stridea, const size_t strideb)
{
/* sample sets of doubles */
size_t i;
const size_t na = 14, nb = 14;
const double rawa[] =
{.0421, .0941, .1064, .0242, .1331,
.0773, .0243, .0815, .1186, .0356,
.0728, .0999, .0614, .0479};
const double rawb[] =
{.1081, .0986, .1566, .1961, .1125,
.1942, .1079, .1021, .1583, .1673,
.1675, .1856, .1688, .1512};
const double raww[] =
{.0000, .0000, .0000, 3.000, .0000,
1.000, 1.000, 1.000, 0.000, .5000,
7.000, 5.000, 4.000, 0.123};
BASE * sorted ;
BASE * groupa = (BASE *) malloc (stridea * na * sizeof(BASE));
BASE * groupb = (BASE *) malloc (strideb * nb * sizeof(BASE));
BASE * w = (BASE *) malloc (strideb * na * sizeof(BASE));
#ifdef BASE_FLOAT
double rel = 1e-6;
#else
double rel = 1e-10;
#endif
for (i = 0 ; i < na ; i++)
groupa[i * stridea] = (BASE) rawa[i] ;
for (i = 0 ; i < na ; i++)
w[i * strideb] = (BASE) raww[i] ;
for (i = 0 ; i < nb ; i++)
groupb[i * strideb] = (BASE) rawb[i] ;
{
double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
double expected = 0.0728;
gsl_test_rel (mean, expected, rel, NAME(gsl_stats) "_mean");
}
{
double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
double var = FUNCTION(gsl_stats,variance_with_fixed_mean) (groupa, stridea, na, mean);
double expected = 0.00113837428571429;
gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_variance_with_fixed_mean");
}
{
double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
double var = FUNCTION(gsl_stats,sd_with_fixed_mean) (groupa, stridea, na, mean);
double expected = 0.0337398026922845;
gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_sd_with_fixed_mean");
}
{
double var = FUNCTION(gsl_stats,variance) (groupb, strideb, nb);
double expected = 0.00124956615384615;
gsl_test_rel (var, expected, rel, NAME(gsl_stats) "_variance");
}
{
double sd = FUNCTION(gsl_stats,sd) (groupa, stridea, na);
double expected = 0.0350134479659107;
gsl_test_rel (sd, expected, rel, NAME(gsl_stats) "_sd");
}
{
double ss = FUNCTION(gsl_stats,tss) (groupb, strideb, nb);
double expected = 0.01624436;
gsl_test_rel (ss, expected, rel, NAME(gsl_stats) "_ss");
}
{
double mean = FUNCTION(gsl_stats,mean) (groupa, stridea, na);
double ss = FUNCTION(gsl_stats,tss_m) (groupa, stridea, na, mean);
double expected = 1.59372400000000e-02;
gsl_test_rel (ss, expected, rel, NAME(gsl_stats) "_ss_m");
}
{
double absdev = FUNCTION(gsl_stats,absdev) (groupa, stridea, na);
double expected = 0.0287571428571429;
gsl_test_rel (absdev, expected, rel, NAME(gsl_stats) "_absdev");
}
{
double skew = FUNCTION(gsl_stats,skew) (groupa, stridea, na);
double expected = 0.0954642051479004;
gsl_test_rel (skew, expected, rel, NAME(gsl_stats) "_skew");
}
{
double kurt = FUNCTION(gsl_stats,kurtosis) (groupa, stridea, na);
double expected = -1.38583851548909 ;
gsl_test_rel (kurt, expected, rel, NAME(gsl_stats) "_kurtosis");
}
{
double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
double expected = 0.0678111523670601;
gsl_test_rel (wmean, expected, rel, NAME(gsl_stats) "_wmean");
}
{
double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
double wvar = FUNCTION(gsl_stats,wvariance_with_fixed_mean) (w, strideb, groupa, stridea, na, wmean);
double expected = 0.000615793060878654;
gsl_test_rel (wvar, expected, rel, NAME(gsl_stats) "_wvariance_with_fixed_mean");
}
{
double est_wvar = FUNCTION(gsl_stats,wvariance) (w, strideb, groupa, stridea, na);
double expected = 0.000769562962860317;
gsl_test_rel (est_wvar, expected, rel, NAME(gsl_stats) "_wvariance");
}
{
double wsd = FUNCTION(gsl_stats,wsd) (w, strideb, groupa, stridea, na);
double expected = 0.0277409978706664;
gsl_test_rel (wsd, expected, rel, NAME(gsl_stats) "_wsd");
}
{
double wtss = FUNCTION(gsl_stats,wtss) (w, strideb, groupa, stridea, na);
double expected = 1.39310864162578e-02;
gsl_test_rel (wtss, expected, rel, NAME(gsl_stats) "_wtss");
}
{
double wmean = FUNCTION(gsl_stats,wmean) (w, strideb, groupa, stridea, na);
double wtss = FUNCTION(gsl_stats,wtss_m) (w, strideb, groupa, stridea, na, wmean);
double expected = 1.39310864162578e-02;
gsl_test_rel (wtss, expected, rel, NAME(gsl_stats) "_wtss_m");
}
{
double wabsdev = FUNCTION(gsl_stats,wabsdev) (w, strideb, groupa, stridea, na);
double expected = 0.0193205027504008;
gsl_test_rel (wabsdev, expected, rel, NAME(gsl_stats) "_wabsdev");
}
{
double wskew = FUNCTION(gsl_stats,wskew) (w, strideb, groupa, stridea, na);
double expected = -0.373631000307076;
gsl_test_rel (wskew, expected, rel, NAME(gsl_stats) "_wskew");
}
{
double wkurt = FUNCTION(gsl_stats,wkurtosis) (w, strideb, groupa, stridea, na);
double expected = -1.48114233353963;
gsl_test_rel (wkurt, expected, rel, NAME(gsl_stats) "_wkurtosis");
}
{
double c = FUNCTION(gsl_stats,covariance) (groupa, stridea, groupb, strideb, nb);
double expected = -0.000139021538461539;
gsl_test_rel (c, expected, rel, NAME(gsl_stats) "_covariance");
}
{
double r = FUNCTION(gsl_stats,correlation) (groupa, stridea, groupb, strideb, nb);
double expected = -0.112322712666074171;
gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_correlation");
}
{
double *work = malloc(2 * na * sizeof(double));
double r = FUNCTION(gsl_stats,spearman) (groupa, stridea, groupb, strideb, nb, work);
double expected = -0.1604395604395604396;
gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_spearman");
free(work);
}
{
double pv = FUNCTION(gsl_stats,pvariance) (groupa, stridea, na, groupb, strideb, nb);
double expected = 0.00123775384615385;
gsl_test_rel (pv, expected, rel, NAME(gsl_stats) "_pvariance");
}
{
double t = FUNCTION(gsl_stats,ttest) (groupa, stridea, na, groupb, strideb, nb);
double expected = -5.67026326985851;
gsl_test_rel (t, expected, rel, NAME(gsl_stats) "_ttest");
}
{
BASE expected = (BASE)0.1331;
gsl_test (FUNCTION(gsl_stats,max) (groupa, stridea, na) != expected,
NAME(gsl_stats) "_max (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
FUNCTION(gsl_stats,max) (groupa, stridea, na), expected);
}
{
BASE min = FUNCTION(gsl_stats,min) (groupa, stridea, na);
BASE expected = (BASE)0.0242;
gsl_test (min != expected,
NAME(gsl_stats) "_min (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
min, expected);
}
{
BASE min, max;
BASE expected_max = (BASE)0.1331;
BASE expected_min = (BASE)0.0242;
FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na);
gsl_test (max != expected_max,
NAME(gsl_stats) "_minmax max (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
max, expected_max);
gsl_test (min != expected_min,
NAME(gsl_stats) "_minmax min (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
min, expected_min);
}
{
int max_index = FUNCTION(gsl_stats,max_index) (groupa, stridea, na);
int expected = 4;
gsl_test (max_index != expected,
NAME(gsl_stats) "_max_index (%d observed vs %d expected)",
max_index, expected);
}
{
int min_index = FUNCTION(gsl_stats,min_index) (groupa, stridea, na);
int expected = 3;
gsl_test (min_index != expected,
NAME(gsl_stats) "_min_index (%d observed vs %d expected)",
min_index, expected);
}
{
size_t min_index, max_index;
size_t expected_max_index = 4;
size_t expected_min_index = 3;
FUNCTION(gsl_stats,minmax_index) (&min_index, &max_index, groupa, stridea, na);
gsl_test (max_index != expected_max_index,
NAME(gsl_stats) "_minmax_index max (%u observed vs %u expected)",
max_index, expected_max_index);
gsl_test (min_index != expected_min_index,
NAME(gsl_stats) "_minmax_index min (%u observed vs %u expected)",
min_index, expected_min_index);
}
sorted = (BASE *) malloc(stridea * na * sizeof(BASE)) ;
for (i = 0 ; i < na ; i++)
sorted[stridea * i] = groupa[stridea * i] ;
TYPE(gsl_sort)(sorted, stridea, na) ;
{
double median = FUNCTION(gsl_stats,median_from_sorted_data)(sorted, stridea, na) ;
double expected = 0.07505;
gsl_test_rel (median,expected, rel,
NAME(gsl_stats) "_median_from_sorted_data (even)");
}
{
double median = FUNCTION(gsl_stats,median_from_sorted_data)(sorted, stridea, na - 1) ;
double expected = 0.0728;
gsl_test_rel (median,expected, rel,
NAME(gsl_stats) "_median_from_sorted_data (odd)");
}
{
BASE * work = (BASE *) malloc (stridea * na * sizeof(BASE));
double expected = 0.07505;
double median;
for (i = 0; i < na; i++)
work[i * stridea] = (BASE) rawa[i];
median = FUNCTION(gsl_stats,median)(work, stridea, na) ;
gsl_test_rel (median, expected, rel,
NAME(gsl_stats) "_median (even)");
free(work);
}
{
BASE * work = (BASE *) malloc (stridea * (na - 1) * sizeof(BASE));
double expected = 0.0728;
double median;
for (i = 0; i < na - 1; i++)
work[i * stridea] = (BASE) sorted[i * stridea];
median = FUNCTION(gsl_stats,median)(work, stridea, na - 1) ;
gsl_test_rel (median, expected, rel,
NAME(gsl_stats) "_median (odd)");
free(work);
}
{
const double trim = 0.2;
double trmean = FUNCTION(gsl_stats,trmean_from_sorted_data)(trim, sorted, stridea, na) ;
double expected = 0.0719;
gsl_test_rel (trmean, expected, rel,
NAME(gsl_stats) "_trmean_from_sorted_data (even)");
}
{
const double trim = 0.2;
double trmean = FUNCTION(gsl_stats,trmean_from_sorted_data)(trim, sorted, stridea, na - 1) ;
double expected = 0.06806666666666666;
gsl_test_rel (trmean, expected, rel,
NAME(gsl_stats) "_trmean_from_sorted_data (odd)");
}
{
double gastwirth = FUNCTION(gsl_stats,gastwirth_from_sorted_data)(sorted, stridea, na) ;
double expected = 0.07271;
gsl_test_rel (gastwirth, expected, rel,
NAME(gsl_stats) "_gastwirth_from_sorted_data (even)");
}
{
double gastwirth = FUNCTION(gsl_stats,gastwirth_from_sorted_data)(sorted, stridea, na - 1) ;
double expected = 0.06794;
gsl_test_rel (gastwirth, expected, rel,
NAME(gsl_stats) "_gastwirth_from_sorted_data (odd)");
}
{
double zeroth = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 0.0) ;
double expected = 0.0242;
gsl_test_rel (zeroth,expected, rel,
NAME(gsl_stats) "_quantile_from_sorted_data (0)");
}
{
double top = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 1.0) ;
double expected = 0.1331;
gsl_test_rel (top,expected, rel,
NAME(gsl_stats) "_quantile_from_sorted_data (100)");
}
{
double median = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na, 0.5) ;
double expected = 0.07505;
gsl_test_rel (median,expected, rel,
NAME(gsl_stats) "_quantile_from_sorted_data (50even)");
}
{
double median = FUNCTION(gsl_stats,quantile_from_sorted_data)(sorted, stridea, na - 1, 0.5);
double expected = 0.0728;
gsl_test_rel (median,expected, rel,
NAME(gsl_stats) "_quantile_from_sorted_data (50odd)");
}
{
size_t k;
BASE * work = (BASE *) malloc (stridea * na * sizeof(BASE));
for (k = 0; k < na; ++k)
{
double expected = sorted[k * stridea];
double kselect;
/* copy rawa[] for each k, since gsl_stats_select() changes input array */
for (i = 0; i < na; i++)
work[i * stridea] = (BASE) rawa[i];
kselect = FUNCTION(gsl_stats,select)(work, stridea, na, k);
gsl_test_rel(kselect, expected, rel, NAME(gsl_stats) "_select");
}
free(work);
}
/* Test for IEEE handling - set third element to NaN */
groupa [3*stridea] = GSL_NAN;
{
BASE max = FUNCTION(gsl_stats,max) (groupa, stridea, na);
BASE expected = GSL_NAN;
gsl_test (!isnan(max),
NAME(gsl_stats) "_max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
max, expected);
}
{
BASE min = FUNCTION(gsl_stats,min) (groupa, stridea, na);
BASE expected = GSL_NAN;
gsl_test (!isnan(min),
NAME(gsl_stats) "_min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
min, expected);
}
{
BASE min, max;
BASE expected_max = GSL_NAN;
BASE expected_min = GSL_NAN;
FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na);
gsl_test (!isnan(max),
NAME(gsl_stats) "_minmax max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
max, expected_max);
gsl_test (!isnan(min),
NAME(gsl_stats) "_minmax min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
min, expected_min);
}
#ifdef FAST
{
BASE min, max;
BASE expected_max = GSL_NAN;
BASE expected_min = GSL_NAN;
FUNCTION(gsl_stats,minmax) (&min, &max, groupa, stridea, na-1);
gsl_test (!isnan(max),
NAME(gsl_stats) "_minmax(-1) max NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
max, expected_max);
gsl_test (!isnan(min),
NAME(gsl_stats) "_minmax(-1) min NaN (" OUT_FORMAT " observed vs " OUT_FORMAT " expected)",
min, expected_min);
}
#endif
{
int max_index = FUNCTION(gsl_stats,max_index) (groupa, stridea, na);
int expected = 3;
gsl_test (max_index != expected,
NAME(gsl_stats) "_max_index NaN (%d observed vs %d expected)",
max_index, expected);
}
{
int min_index = FUNCTION(gsl_stats,min_index) (groupa, stridea, na);
int expected = 3;
gsl_test (min_index != expected,
NAME(gsl_stats) "_min_index NaN (%d observed vs %d expected)",
min_index, expected);
}
{
size_t min_index, max_index;
size_t expected_max_index = 3;
size_t expected_min_index = 3;
FUNCTION(gsl_stats,minmax_index) (&min_index, &max_index, groupa, stridea, na);
gsl_test (max_index != expected_max_index,
NAME(gsl_stats) "_minmax_index max NaN (%u observed vs %u expected)",
max_index, expected_max_index);
gsl_test (min_index != expected_min_index,
NAME(gsl_stats) "_minmax_index min NaN (%u observed vs %u expected)",
min_index, expected_min_index);
}
/* restore */
groupa [3*stridea] = (BASE) rawa[3];
{
double * work = (double *) malloc (na * sizeof(double));
double expected = 0.02925;
double mad0 = FUNCTION(gsl_stats,mad0)(groupa, stridea, na, work);
gsl_test_rel (mad0, expected, rel, NAME(gsl_stats) "_mad0 (even)");
free(work);
}
{
double * work = (double *) malloc ((na - 1) * sizeof(double));
double expected = 0.02910;
double mad0 = FUNCTION(gsl_stats,mad0)(groupa, stridea, na - 1, work);
gsl_test_rel (mad0, expected, rel, NAME(gsl_stats) "_mad0 (odd)");
free(work);
}
{
BASE *work = malloc(na * sizeof(BASE));
double r = FUNCTION(gsl_stats,Sn_from_sorted_data) (sorted, stridea, na, work);
double expected = 0.04007136;
gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Sn_from_sorted_data (even)");
free(work);
}
{
BASE *work = malloc((na - 1) * sizeof(BASE));
double r = FUNCTION(gsl_stats,Sn_from_sorted_data) (sorted, stridea, na - 1, work);
double expected = 0.03728599834710744;
gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Sn_from_sorted_data (odd)");
free(work);
}
{
BASE *work = malloc(3 * na * sizeof(BASE));
int *work_int = malloc(5 * na * sizeof(int));
double r = FUNCTION(gsl_stats,Qn_from_sorted_data) (sorted, stridea, na, work, work_int);
double expected = 0.04113672759664409;
gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Qn_from_sorted_data (even)");
free(work);
free(work_int);
}
{
BASE *work = malloc(3 * (na - 1) * sizeof(BASE));
int *work_int = malloc(5 * (na - 1) * sizeof(int));
double r = FUNCTION(gsl_stats,Qn_from_sorted_data) (sorted, stridea, na - 1, work, work_int);
double expected = 0.03684305546303433;
gsl_test_rel (r, expected, rel, NAME(gsl_stats) "_Qn_from_sorted_data (odd)");
free(work);
free(work_int);
}
free (sorted);
free (groupa);
free (groupb);
free (w);
}