diff options
Diffstat (limited to 'stats.h')
| -rw-r--r-- | stats.h | 213 |
1 files changed, 213 insertions, 0 deletions
@@ -0,0 +1,213 @@ +#include "c23-fallback.h" +#include <tgmath.h> +#include <stdio.h> + +/** + ** @file + ** @brief Collect some simple statistics online as we go. + ** + ** This uses a generalization of Welford's trick to compute running + ** mean and variance. See + ** + ** Philippe Pébay. Formulas for robust, one-pass parallel computation + ** of covariances and arbitrary-order statistical moments. Technical + ** Report SAND2008-6212, SANDIA, 2008. URL + ** http://prod.sandia.gov/techlib/access-control.cgi/2008/086212.pdf. + **/ + +/** + ** @brief A simple data structure to collect the 0th to 3rd moment of + ** a statistic. + ** + ** @warning Since this also uses a @c double for the number of + ** samples, the validity of all this is restricted to about + ** @f$2^{50} \approx 10^{15}@f$ samples. + **/ +struct stats { + double moment[4]; +}; +typedef struct stats stats; + + +/** + ** @brief Return the number of samples that had been entered into the + ** statistic @a c. + **/ +inline +double stats_samples(stats c[static 1]) [[__unsequenced__]] { + return c->moment[0]; +} + +/** + ** @brief Return the mean value of the samples that had been entered + ** into the statistic @a c. + **/ +inline +double stats_mean(stats c[static 1]) [[__unsequenced__]] { + return c->moment[1]; +} + +/** + ** @brief Return the variance of the samples that had been entered + ** into the statistic @a c. + **/ +inline +double stats_var(stats c[static 1]) [[__reproducible__]] { + return c->moment[2]/stats_samples(c); +} + +// Clang has this notorious design flaw in <tgmath.h> where they +// declare interfaces as being static. In some cases this results in +// an avalanche of confusing messages of static functions that would +// be used in an inline function. +#if __clang_major__ +#pragma clang diagnostic push +#pragma clang diagnostic ignored "-Wstatic-in-inline" +#endif + +/** + ** @brief Return the standard deviation of the samples that had been + ** entered into the statistic @a c. + **/ +inline +double stats_sdev(stats c[static 1]) [[__reproducible__]] { + return sqrt(stats_var(c)); +} + +/** + ** @brief Return the relative standard deviation of the samples that + ** had been entered into the statistic @a c. + **/ +inline +double stats_rsdev(stats c[static 1]) [[__reproducible__]] { + return sqrt(stats_var(c))/stats_mean(c); +} + +/** + ** @brief Return the normalized skew of the samples that had been + ** entered into the statistic @a c. + **/ +inline +double stats_skew(stats c[static 1]) [[__reproducible__]] { + double var = stats_var(c); + return (c->moment[3]/pow(var, 1.5))/stats_samples(c); +} + +/** + ** @brief Return the unbiased variance of the samples that had been + ** entered into the statistic @a c. + ** + ** Use Bessel's correction to have an estimation of the unbiased + ** variance of the overall population. + **/ +inline +double stats_var_unbiased(stats c[static 1]) [[__reproducible__]] { + return c->moment[2]/(stats_samples(c)-1); +} + +/** + ** @brief Return the unbiased standard deviation of the samples that + ** had been entered into the statistic @a c. + ** + ** Use Bessel's correction to have an less biased estimation of the + ** variance of the overall population. + **/ +inline +double stats_sdev_unbiased(stats c[static 1]) [[__reproducible__]] { + return sqrt(stats_var_unbiased(c)); +} + +#if __clang_major__ +#pragma clang diagnostic pop +#endif + + + +/** + ** @brief Return the unbiased relative standard deviation of the + ** samples that had been entered into the statistic @a c. + **/ +inline +double stats_rsdev_unbiased(stats c[static 1]) [[__reproducible__]] { + return stats_rsdev(c)*(1+1/(4*stats_samples(c))); +} + +/** + ** @brief Add value @a val to the statistic @a c. + ** + ** @c moments is the number of statistic moments that is collected, + ** it has to be between `0` and `3`, including. + **/ +inline +void stats_collect(stats c[static 1], + double val, unsigned moments) + [[__reproducible__]] { + double n = stats_samples(c); + double n0 = n-1; + double n1 = n+1; + double delta0 = 1; + double delta = val - stats_mean(c); + double delta1 = delta/n1; + double delta2 = delta1*delta*n; + switch (moments) { + default: + c->moment[3] += (delta2*n0 - 3*c->moment[2])*delta1; + [[__fallthrough__]]; + case 2: + c->moment[2] += delta2; + [[__fallthrough__]]; + case 1: + c->moment[1] += delta1; + [[__fallthrough__]]; + case 0: + c->moment[0] += delta0; + } +} + +/** + ** @brief Add value @a val to the statistic @a c. + ** + ** Only the number of samples is collected. + **/ +inline +void stats_collect0(stats c[static 1], + double val) + [[__reproducible__]] { + stats_collect(c, val, 0); +} + +/** + ** @brief Add value @a val to the statistic @a c. + ** + ** Only the number of samples and mean is collected. + **/ +inline +void stats_collect1(stats c[static 1], + double val) + [[__reproducible__]] { + stats_collect(c, val, 1); +} + +/** + ** @brief Add value @a val to the statistic @a c. + ** + ** The number of samples, mean and standard deviation is collected. + **/ +inline +void stats_collect2(stats c[static 1], + double val) + [[__reproducible__]] { + stats_collect(c, val, 2); +} + +/** + ** @brief Add value @a val to the statistic @a c. + ** + ** All statistics up to the skew are collected. + **/ +inline +void stats_collect3(stats c[static 1], + double val) + [[__reproducible__]] { + stats_collect(c, val, 3); +} |