diff options
Diffstat (limited to 'src/statistical.c')
| -rw-r--r-- | src/statistical.c | 142 |
1 files changed, 71 insertions, 71 deletions
diff --git a/src/statistical.c b/src/statistical.c index 807bc25ea..d4b33fd5a 100644 --- a/src/statistical.c +++ b/src/statistical.c @@ -2,7 +2,7 @@ // -------------------------------------------------------------------------------------------------------------------- -inline long double sum_and_count(long double *series, size_t entries, size_t *count) { +inline LONG_DOUBLE sum_and_count(const LONG_DOUBLE *series, size_t entries, size_t *count) { if(unlikely(entries == 0)) { if(likely(count)) *count = 0; @@ -18,10 +18,10 @@ inline long double sum_and_count(long double *series, size_t entries, size_t *co } size_t i, c = 0; - long double sum = 0; + LONG_DOUBLE sum = 0; for(i = 0; i < entries ; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; c++; sum += value; @@ -36,44 +36,44 @@ inline long double sum_and_count(long double *series, size_t entries, size_t *co return sum; } -inline long double sum(long double *series, size_t entries) { +inline LONG_DOUBLE sum(const LONG_DOUBLE *series, size_t entries) { return sum_and_count(series, entries, NULL); } -inline long double average(long double *series, size_t entries) { +inline LONG_DOUBLE average(const LONG_DOUBLE *series, size_t entries) { size_t count = 0; - long double sum = sum_and_count(series, entries, &count); + LONG_DOUBLE sum = sum_and_count(series, entries, &count); if(unlikely(count == 0)) return NAN; - return sum / (long double)count; + return sum / (LONG_DOUBLE)count; } // -------------------------------------------------------------------------------------------------------------------- -long double moving_average(long double *series, size_t entries, size_t period) { +LONG_DOUBLE moving_average(const LONG_DOUBLE *series, size_t entries, size_t period) { if(unlikely(period <= 0)) return 0.0; size_t i, count; - long double sum = 0, avg = 0; - long double p[period]; + LONG_DOUBLE sum = 0, avg = 0; + LONG_DOUBLE p[period]; for(count = 0; count < period ; count++) p[count] = 0.0; for(i = 0, count = 0; i < entries; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; if(unlikely(count < period)) { sum += value; - avg = (count == period - 1) ? sum / (long double)period : 0; + avg = (count == period - 1) ? sum / (LONG_DOUBLE)period : 0; } else { sum = sum - p[count % period] + value; - avg = sum / (long double)period; + avg = sum / (LONG_DOUBLE)period; } p[count % period] = value; @@ -86,8 +86,8 @@ long double moving_average(long double *series, size_t entries, size_t period) { // -------------------------------------------------------------------------------------------------------------------- static int qsort_compare(const void *a, const void *b) { - long double *p1 = (long double *)a, *p2 = (long double *)b; - long double n1 = *p1, n2 = *p2; + LONG_DOUBLE *p1 = (LONG_DOUBLE *)a, *p2 = (LONG_DOUBLE *)b; + LONG_DOUBLE n1 = *p1, n2 = *p2; if(unlikely(isnan(n1) || isnan(n2))) { if(isnan(n1) && !isnan(n2)) return -1; @@ -105,17 +105,17 @@ static int qsort_compare(const void *a, const void *b) { return 0; } -inline void sort_series(long double *series, size_t entries) { - qsort(series, entries, sizeof(long double), qsort_compare); +inline void sort_series(LONG_DOUBLE *series, size_t entries) { + qsort(series, entries, sizeof(LONG_DOUBLE), qsort_compare); } -inline long double *copy_series(long double *series, size_t entries) { - long double *copy = mallocz(sizeof(long double) * entries); - memcpy(copy, series, sizeof(long double) * entries); +inline LONG_DOUBLE *copy_series(const LONG_DOUBLE *series, size_t entries) { + LONG_DOUBLE *copy = mallocz(sizeof(LONG_DOUBLE) * entries); + memcpy(copy, series, sizeof(LONG_DOUBLE) * entries); return copy; } -long double median_on_sorted_series(long double *series, size_t entries) { +LONG_DOUBLE median_on_sorted_series(const LONG_DOUBLE *series, size_t entries) { if(unlikely(entries == 0)) return NAN; @@ -125,7 +125,7 @@ long double median_on_sorted_series(long double *series, size_t entries) { if(unlikely(entries == 2)) return (series[0] + series[1]) / 2; - long double avg; + LONG_DOUBLE avg; if(entries % 2 == 0) { size_t m = entries / 2; avg = (series[m] + series[m + 1]) / 2; @@ -137,7 +137,7 @@ long double median_on_sorted_series(long double *series, size_t entries) { return avg; } -long double median(long double *series, size_t entries) { +LONG_DOUBLE median(const LONG_DOUBLE *series, size_t entries) { if(unlikely(entries == 0)) return NAN; @@ -147,10 +147,10 @@ long double median(long double *series, size_t entries) { if(unlikely(entries == 2)) return (series[0] + series[1]) / 2; - long double *copy = copy_series(series, entries); + LONG_DOUBLE *copy = copy_series(series, entries); sort_series(copy, entries); - long double avg = median_on_sorted_series(copy, entries); + LONG_DOUBLE avg = median_on_sorted_series(copy, entries); freez(copy); return avg; @@ -158,18 +158,18 @@ long double median(long double *series, size_t entries) { // -------------------------------------------------------------------------------------------------------------------- -long double moving_median(long double *series, size_t entries, size_t period) { +LONG_DOUBLE moving_median(const LONG_DOUBLE *series, size_t entries, size_t period) { if(entries <= period) return median(series, entries); - long double *data = copy_series(series, entries); + LONG_DOUBLE *data = copy_series(series, entries); size_t i; for(i = period; i < entries; i++) { data[i - period] = median(&series[i - period], period); } - long double avg = median(data, entries - period); + LONG_DOUBLE avg = median(data, entries - period); freez(data); return avg; } @@ -177,13 +177,13 @@ long double moving_median(long double *series, size_t entries, size_t period) { // -------------------------------------------------------------------------------------------------------------------- // http://stackoverflow.com/a/15150143/4525767 -long double running_median_estimate(long double *series, size_t entries) { - long double median = 0.0f; - long double average = 0.0f; +LONG_DOUBLE running_median_estimate(const LONG_DOUBLE *series, size_t entries) { + LONG_DOUBLE median = 0.0f; + LONG_DOUBLE average = 0.0f; size_t i; for(i = 0; i < entries ; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; average += ( value - average ) * 0.1f; // rough running average. @@ -195,7 +195,7 @@ long double running_median_estimate(long double *series, size_t entries) { // -------------------------------------------------------------------------------------------------------------------- -long double standard_deviation(long double *series, size_t entries) { +LONG_DOUBLE standard_deviation(const LONG_DOUBLE *series, size_t entries) { if(unlikely(entries < 1)) return NAN; @@ -203,10 +203,10 @@ long double standard_deviation(long double *series, size_t entries) { return series[0]; size_t i, count = 0; - long double sum = 0; + LONG_DOUBLE sum = 0; for(i = 0; i < entries ; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; count++; @@ -219,10 +219,10 @@ long double standard_deviation(long double *series, size_t entries) { if(unlikely(count == 1)) return sum; - long double average = sum / (long double)count; + LONG_DOUBLE average = sum / (LONG_DOUBLE)count; for(i = 0, count = 0, sum = 0; i < entries ; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; count++; @@ -235,29 +235,29 @@ long double standard_deviation(long double *series, size_t entries) { if(unlikely(count == 1)) return average; - long double variance = sum / (long double)(count - 1); // remove -1 to have a population stddev + LONG_DOUBLE variance = sum / (LONG_DOUBLE)(count - 1); // remove -1 to have a population stddev - long double stddev = sqrtl(variance); + LONG_DOUBLE stddev = sqrtl(variance); return stddev; } // -------------------------------------------------------------------------------------------------------------------- -long double single_exponential_smoothing(long double *series, size_t entries, long double alpha) { +LONG_DOUBLE single_exponential_smoothing(const LONG_DOUBLE *series, size_t entries, LONG_DOUBLE alpha) { size_t i, count = 0; - long double level = 0, sum = 0; + LONG_DOUBLE level = 0, sum = 0; if(unlikely(isnan(alpha))) alpha = 0.3; for(i = 0; i < entries ; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; count++; sum += value; - long double last_level = level; + LONG_DOUBLE last_level = level; level = alpha * value + (1.0 - alpha) * last_level; } @@ -267,9 +267,9 @@ long double single_exponential_smoothing(long double *series, size_t entries, lo // -------------------------------------------------------------------------------------------------------------------- // http://grisha.org/blog/2016/02/16/triple-exponential-smoothing-forecasting-part-ii/ -long double double_exponential_smoothing(long double *series, size_t entries, long double alpha, long double beta, long double *forecast) { +LONG_DOUBLE double_exponential_smoothing(const LONG_DOUBLE *series, size_t entries, LONG_DOUBLE alpha, LONG_DOUBLE beta, LONG_DOUBLE *forecast) { size_t i, count = 0; - long double level = series[0], trend, sum; + LONG_DOUBLE level = series[0], trend, sum; if(unlikely(isnan(alpha))) alpha = 0.3; @@ -285,13 +285,13 @@ long double double_exponential_smoothing(long double *series, size_t entries, lo sum = series[0]; for(i = 1; i < entries ; i++) { - long double value = series[i]; + LONG_DOUBLE value = series[i]; if(unlikely(isnan(value) || isinf(value))) continue; count++; sum += value; - long double last_level = level; + LONG_DOUBLE last_level = level; level = alpha * value + (1.0 - alpha) * (level + trend); trend = beta * (level - last_level) + (1.0 - beta) * trend; @@ -327,24 +327,24 @@ long double double_exponential_smoothing(long double *series, size_t entries, lo * s[t] = γ (Y[t] / a[t]) + (1-γ) s[t-p] */ static int __HoltWinters( - long double *series, + const LONG_DOUBLE *series, int entries, // start_time + h - long double alpha, // alpha parameter of Holt-Winters Filter. - long double beta, // beta parameter of Holt-Winters Filter. If set to 0, the function will do exponential smoothing. - long double gamma, // gamma parameter used for the seasonal component. If set to 0, an non-seasonal model is fitted. + LONG_DOUBLE alpha, // alpha parameter of Holt-Winters Filter. + LONG_DOUBLE beta, // beta parameter of Holt-Winters Filter. If set to 0, the function will do exponential smoothing. + LONG_DOUBLE gamma, // gamma parameter used for the seasonal component. If set to 0, an non-seasonal model is fitted. - int *seasonal, - int *period, - long double *a, // Start value for level (a[0]). - long double *b, // Start value for trend (b[0]). - long double *s, // Vector of start values for the seasonal component (s_1[0] ... s_p[0]) + const int *seasonal, + const int *period, + const LONG_DOUBLE *a, // Start value for level (a[0]). + const LONG_DOUBLE *b, // Start value for trend (b[0]). + LONG_DOUBLE *s, // Vector of start values for the seasonal component (s_1[0] ... s_p[0]) /* return values */ - long double *SSE, // The final sum of squared errors achieved in optimizing - long double *level, // Estimated values for the level component (size entries - t + 2) - long double *trend, // Estimated values for the trend component (size entries - t + 2) - long double *season // Estimated values for the seasonal component (size entries - t + 2) + LONG_DOUBLE *SSE, // The final sum of squared errors achieved in optimizing + LONG_DOUBLE *level, // Estimated values for the level component (size entries - t + 2) + LONG_DOUBLE *trend, // Estimated values for the trend component (size entries - t + 2) + LONG_DOUBLE *season // Estimated values for the seasonal component (size entries - t + 2) ) { if(unlikely(entries < 4)) @@ -352,13 +352,13 @@ static int __HoltWinters( int start_time = 2; - long double res = 0, xhat = 0, stmp = 0; + LONG_DOUBLE res = 0, xhat = 0, stmp = 0; int i, i0, s0; /* copy start values to the beginning of the vectors */ level[0] = *a; if(beta > 0) trend[0] = *b; - if(gamma > 0) memcpy(season, s, *period * sizeof(long double)); + if(gamma > 0) memcpy(season, s, *period * sizeof(LONG_DOUBLE)); for(i = start_time - 1; i < entries; i++) { /* indices for period i */ @@ -404,7 +404,7 @@ static int __HoltWinters( return 1; } -long double holtwinters(long double *series, size_t entries, long double alpha, long double beta, long double gamma, long double *forecast) { +LONG_DOUBLE holtwinters(const LONG_DOUBLE *series, size_t entries, LONG_DOUBLE alpha, LONG_DOUBLE beta, LONG_DOUBLE gamma, LONG_DOUBLE *forecast) { if(unlikely(isnan(alpha))) alpha = 0.3; @@ -416,15 +416,15 @@ long double holtwinters(long double *series, size_t entries, long double alpha, int seasonal = 0; int period = 0; - long double a0 = series[0]; - long double b0 = 0; - long double s[] = {}; + LONG_DOUBLE a0 = series[0]; + LONG_DOUBLE b0 = 0; + LONG_DOUBLE s[] = {}; - long double errors = 0.0; + LONG_DOUBLE errors = 0.0; size_t nb_computations = entries; - long double *estimated_level = callocz(nb_computations, sizeof(long double)); - long double *estimated_trend = callocz(nb_computations, sizeof(long double)); - long double *estimated_season = callocz(nb_computations, sizeof(long double)); + LONG_DOUBLE *estimated_level = callocz(nb_computations, sizeof(LONG_DOUBLE)); + LONG_DOUBLE *estimated_trend = callocz(nb_computations, sizeof(LONG_DOUBLE)); + LONG_DOUBLE *estimated_season = callocz(nb_computations, sizeof(LONG_DOUBLE)); int ret = __HoltWinters( series, @@ -443,7 +443,7 @@ long double holtwinters(long double *series, size_t entries, long double alpha, estimated_season ); - long double value = estimated_level[nb_computations - 1]; + LONG_DOUBLE value = estimated_level[nb_computations - 1]; if(forecast) *forecast = 0.0; |