diff options
Diffstat (limited to '')
| -rw-r--r-- | test/unit/regress.c | 119 |
1 files changed, 119 insertions, 0 deletions
diff --git a/test/unit/regress.c b/test/unit/regress.c new file mode 100644 index 0000000..f47d1c4 --- /dev/null +++ b/test/unit/regress.c @@ -0,0 +1,119 @@ +/* + ********************************************************************** + * Copyright (C) Miroslav Lichvar 2017 + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of version 2 of the GNU General Public License as + * published by the Free Software Foundation. + * + * 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. + * + ********************************************************************** + */ +#include <regress.c> +#include "test.h" + +#define POINTS 64 + +void +test_unit(void) +{ + double x[POINTS], x2[POINTS], y[POINTS], w[POINTS]; + double b0, b1, b2, s2, sb0, sb1, slope, slope2, intercept, sd, median; + double xrange, yrange, wrange, x2range; + int i, j, n, m, c1, c2, c3, runs, best_start, dof; + + for (n = 3; n <= POINTS; n++) { + for (i = 0; i < 200; i++) { + slope = TST_GetRandomDouble(-0.1, 0.1); + intercept = TST_GetRandomDouble(-1.0, 1.0); + sd = TST_GetRandomDouble(1e-6, 1e-4); + slope2 = (random() % 2 ? 1 : -1) * TST_GetRandomDouble(0.1, 0.5); + + DEBUG_LOG("iteration %d n=%d intercept=%e slope=%e sd=%e", + i, n, intercept, slope, sd); + + for (j = 0; j < n; j++) { + x[j] = -j; + y[j] = intercept + slope * x[j] + (j % 2 ? 1 : -1) * TST_GetRandomDouble(1e-6, sd); + w[j] = TST_GetRandomDouble(1.0, 2.0); + x2[j] = (y[j] - intercept - slope * x[j]) / slope2; + } + + RGR_WeightedRegression(x, y, w, n, &b0, &b1, &s2, &sb0, &sb1); + DEBUG_LOG("WR b0=%e b1=%e s2=%e sb0=%e sb1=%e", b0, b1, s2, sb0, sb1); + TEST_CHECK(fabs(b0 - intercept) < sd + 1e-3); + TEST_CHECK(fabs(b1 - slope) < sd); + + if (RGR_FindBestRegression(x, y, w, n, 0, 3, &b0, &b1, &s2, &sb0, &sb1, + &best_start, &runs, &dof)) { + DEBUG_LOG("BR b0=%e b1=%e s2=%e sb0=%e sb1=%e runs=%d bs=%d dof=%d", + b0, b1, s2, sb0, sb1, runs, best_start, dof); + + TEST_CHECK(fabs(b0 - intercept) < sd + 1e-3); + TEST_CHECK(fabs(b1 - slope) < sd); + } + + if (RGR_MultipleRegress(x, x2, y, n, &b2)) { + DEBUG_LOG("MR b2=%e", b2); + TEST_CHECK(fabs(b2 - slope2) < 1e-6); + } + + for (j = 0; j < n / 7; j++) + y[random() % n] += 100 * sd; + + if (RGR_FindBestRobustRegression(x, y, n, 1e-8, &b0, &b1, &runs, &best_start)) { + DEBUG_LOG("BRR b0=%e b1=%e runs=%d bs=%d", b0, b1, runs, best_start); + + TEST_CHECK(fabs(b0 - intercept) < sd + 1e-2); + TEST_CHECK(fabs(b1 - slope) < 5.0 * sd); + } + + for (j = 0; j < n; j++) + x[j] = random() % 4 * TST_GetRandomDouble(-1000, 1000); + + median = RGR_FindMedian(x, n); + + for (j = c1 = c2 = c3 = 0; j < n; j++) { + if (x[j] < median) + c1++; + if (x[j] > median) + c3++; + else + c2++; + } + + TEST_CHECK(c1 + c2 >= c3 && c1 <= c2 + c3); + + xrange = TST_GetRandomDouble(1e-6, pow(10.0, random() % 10)); + yrange = random() % 3 * TST_GetRandomDouble(0.0, pow(10.0, random() % 10)); + wrange = random() % 3 * TST_GetRandomDouble(0.0, pow(10.0, random() % 10)); + x2range = random() % 3 * TST_GetRandomDouble(0.0, pow(10.0, random() % 10)); + m = random() % n; + + for (j = 0; j < n; j++) { + x[j] = (j ? x[j - 1] : 0.0) + TST_GetRandomDouble(1e-6, xrange); + y[j] = TST_GetRandomDouble(-yrange, yrange); + w[j] = 1.0 + TST_GetRandomDouble(0.0, wrange); + x2[j] = TST_GetRandomDouble(-x2range, x2range); + } + + RGR_WeightedRegression(x, y, w, n, &b0, &b1, &s2, &sb0, &sb1); + + if (RGR_FindBestRegression(x + m, y + m, w, n - m, m, 3, &b0, &b1, &s2, &sb0, &sb1, + &best_start, &runs, &dof)) + ; + if (RGR_MultipleRegress(x, x2, y, n, &b2)) + ; + if (RGR_FindBestRobustRegression(x, y, n, 1e-8, &b0, &b1, &runs, &best_start)) + ; + } + } +} |