From db12a71fccf0bfc35d8d34ba9bde94b612b4d826 Mon Sep 17 00:00:00 2001 From: borchehq Date: Fri, 16 Aug 2024 10:40:04 +0200 Subject: [PATCH] Fixed a bug in rstats_central_moment_finalize where 1st order moments would always be provided regardless of p. --- src/rstats/include/rstats.h | 26 ++++++++++++++------------ src/rstats/test/test_rstats.c | 32 ++++++++++++++++++++++++++------ 2 files changed, 40 insertions(+), 18 deletions(-) diff --git a/src/rstats/include/rstats.h b/src/rstats/include/rstats.h index 4bc8c34..e786ca8 100644 --- a/src/rstats/include/rstats.h +++ b/src/rstats/include/rstats.h @@ -314,8 +314,8 @@ inline void rstats_central_moment(double x, double w, double *buffer, uint64_t p * * @param results A pointer to an array of doubles where the final mean and * central moments will be stored: - * - `results[0]` will store the final mean value. - * - `results[1]` to `results[p+1]` will store the final central + * - `results[p + 1]` will store the final mean value. + * - `results[0]` to `results[p]` will store the final central * moments from the 0th to the p-th order. * results must point to an array of length p + 2. * @param buffer A pointer to an array of doubles used in @@ -328,20 +328,22 @@ inline void rstats_central_moment(double x, double w, double *buffer, uint64_t p */ inline void rstats_central_moment_finalize(double *results, double *buffer, uint64_t p, bool standardize) { - results[0] = buffer[1]; // Mean. - results[1] = 1.0; // 0th central moment is always 1. - results[2] = 0.0; // 1st central moment is always 0. - for(uint64_t i = 3; i < p + 2; i++) { - results[i] = buffer[i - 1] / buffer[0]; + for(uint64_t i = 2; i < p + 1; i++) { + results[i] = buffer[i] / buffer[0]; } if(standardize) { - results[1] = 1.0; // 0th standardized central moment is always 1. - results[2] = 0.0; // 1st standardized central moment is always 0. - for(uint64_t i = 4; i < p + 2; i++) { - results[i] = results[i] / rstats_pow(sqrt(results[3]), i - 1); + for(uint64_t i = 3; i < p + 1; i++) { + results[i] = results[i] / rstats_pow(sqrt(results[2]), i); + } + if(p >= 2) { + results[2] = 1.0; // 2nd standardized central moment is always 1. } - results[3] = 1.0; // 2nd standardized central moment is always 1. } + results[0] = 1.0; // 0th standardized central moment is always 1. + if(p >= 1) { + results[1] = 0.0; // 1st standardized central moment is always 0. + } + results[p + 1] = buffer[1]; // Mean. } /** diff --git a/src/rstats/test/test_rstats.c b/src/rstats/test/test_rstats.c index 2b2541d..b8e0b29 100644 --- a/src/rstats/test/test_rstats.c +++ b/src/rstats/test/test_rstats.c @@ -117,6 +117,7 @@ void test_central_moment() { double results_2[4] = {0.0}, buffer_2[5] = {0.0}; double sum_weights = 0.0, sum_weights_comp = 0.0; double tmp = 0.0; + double results_3[3], buffer_3[3]; for(size_t i = 0; i < 10; i++) { rstats_central_moment(x[i], weights[i], buffer, p); @@ -126,14 +127,14 @@ void test_central_moment() { mass += weights[i] * x[i]; sum_weights_comp += weights[i]; mean_cmp = mass / sum_weights_comp; - //printf("%f, %f\n", results[0], wmoments_comp[0] / sum_weights_comp); - assert(fabs(results[0] - mean_cmp) < 1e-7); + //printf("%f, %f\n", results[16], mean_cmp); + assert(fabs(results[16] - mean_cmp) < 1e-7); for(size_t k = 0; k < p; k++) { for(size_t j = 0; j < i + 1; j++) { wmoments_comp[k] += weights[j] * rstats_pow(x[j] - mean_cmp, k); } - //printf("%f, %f\n", results[k + 1], wmoments_comp[k] / sum_weights_comp); - assert(fabs(results[k + 1] - wmoments_comp[k] / sum_weights_comp) < 1e-2); + //printf("%f, %f\n", results[k], wmoments_comp[k] / sum_weights_comp); + assert(fabs(results[k] - wmoments_comp[k] / sum_weights_comp) < 1e-2); wmoments_comp[k] = 0.0; } } @@ -148,8 +149,27 @@ void test_central_moment() { rstats_kurtosis(x[i], weights[i], buffer_2); rstats_kurtosis_finalize(results_2, buffer_2); if(i > 0) { - assert(fabs(results[4] - results_2[2]) < 1e-2); - assert(fabs(results[5] - results_2[3]) < 1e-2); + assert(fabs(results[3] - results_2[2]) < 1e-2); + assert(fabs(results[4] - results_2[3]) < 1e-2); + } + } + + for(size_t k = 0; k < 3; k++) { + for(size_t i = 0; i < 10; i++) { + rstats_central_moment(x[i], weights[i], buffer_3, k); + rstats_central_moment_finalize(results_3, buffer_3, k, false); + assert(results_3[0] == 1.0); + if(k >= 1) { + assert(results_3[1] == 0.0); + } + rstats_central_moment_finalize(results_3, buffer_3, k, true); + assert(results_3[0] == 1.0); + if(k >= 1) { + assert(results_3[1] == 0.0); + } + if(k >= 2) { + assert(results_3[2] == 1.0); + } } } }