Added RunStats unit tests comparing statistics to textbook definitions, comparing the underlying data to an independent implementation and checking the combination operation

Fixed a subtle sign error in the combination operation of two RunStats objects that resulted in incorrect values for the skewness and kurtosis of the result

Author: giltirn

include/chimbuko/util/RunStats.hpp

@@ -24,13 +24,14 @@ namespace chimbuko {
     /**
      * @brief Internal state of RunStats object
      *
+     * Note the variables in https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance are M2,M3,M4. The mappings are provided in the comments below.
      */
     struct State {
       double count; /**< count of instances */
       double eta; /**< mean */
-      double rho;
-      double tau;
-      double phi;
+      double rho; /**< = M2 = \sum_i (x_i - \bar x)^2 */
+      double tau; /**< = M3 = \sum_i (x_i - \bar x)^3 */
+      double phi; /**< = M4 = \sum_i (x_i - \bar x)^4 */
       double min; /**< minimum */
       double max; /**< maximum */
       double acc; /**< sum */
@@ -231,6 +232,13 @@ namespace chimbuko {
      */
     double accumulate() const;
     double mean() const;
+
+    /**
+     * @brief Return the variance of the data
+     *
+     * If ddof=1 (default) the variance will include Bessel's correction, and represents an estimate of the population variance.
+     * If ddof=0 the variance will be the variance of the sample
+     */
     double variance(double ddof=1.0) const;
     double stddev(double ddof=1.0) const;
     double skewness() const;
@@ -241,6 +249,11 @@ namespace chimbuko {
      */
     void set_do_accumulate(bool do_accumulate) { m_do_accumulate = do_accumulate; }
 
+    /**
+     * @brief Determine whether the sum of all values is to be maintained
+     */
+    bool get_do_accumulate() const{ return m_do_accumulate; }
+
     /**
      * @brief Get the current statistics as a JSON object
      */
@@ -255,7 +268,7 @@ namespace chimbuko {
     /**
      * @brief Combine two RunStats instances such that the resulting statistics are the union of the two
      */
-    friend RunStats operator+(const RunStats a, const RunStats b);
+    friend RunStats operator+(const RunStats &a, const RunStats &b);
 
     /**
      * @brief Combine two RunStats instances such that the resulting statistics are the union of the two
@@ -277,7 +290,7 @@ namespace chimbuko {
     bool m_do_accumulate; /**< True if the sum of the input values are maintained */
   };
 
-  RunStats operator+(const RunStats a, const RunStats b);
+  RunStats operator+(const RunStats &a, const RunStats &b);
   bool operator==(const RunStats& a, const RunStats& b);
   bool operator!=(const RunStats& a, const RunStats& b);
 

src/util/RunStats.cpp

@@ -155,13 +155,13 @@ double RunStats::kurtosis() const {
     return m_state.count * m_state.phi / (m_state.rho * m_state.rho) - 3.0;
 }
 
-RunStats chimbuko::operator+(const RunStats a, const RunStats b)
+RunStats chimbuko::operator+(const RunStats &a, const RunStats &b)
 {
     double sum_count = a.m_state.count + b.m_state.count;
     if (sum_count == 0.0)
         return RunStats();
 
-    double delta  = a.m_state.eta - b.m_state.eta;
+    double delta  = b.m_state.eta - a.m_state.eta;
     double delta2 = delta * delta;
     double delta3 = delta * delta2;
     double delta4 = delta2 * delta2;

test/unit_tests/util/RunStats.cpp

@@ -2,9 +2,266 @@
 #include "gtest/gtest.h"
 #include <cereal/archives/portable_binary.hpp>
 #include <sstream>
+#include "../unit_test_common.hpp"
 
 using namespace chimbuko;
 
+//Textbook definitions
+double mean(const std::vector<double> &a){
+  double r=0;
+  for(double v:a) r+=v;
+  return r/a.size();
+}
+double variance(const std::vector<double> &a, bool incl_bessel = true){
+  double n = a.size();
+  double r=0,r2=0;
+  for(double v:a){ r+=v; r2+=v*v; }
+  r = (r2/n - r/n*r/n);
+  if(incl_bessel) r *= n/(n-1);  //include Bessel's correction by default
+  return r;
+}
+double skewness(const std::vector<double> &a){
+  double mu = mean(a);
+  double sigma = sqrt(variance(a,false));
+  double r = 0;
+  for(double v:a) r += pow( (v - mu)/sigma, 3 );
+  return r/a.size();
+}
+double kurtosis(const std::vector<double> &a){ //technically excess kurtosis
+  double mu = mean(a);
+  double sigma = sqrt(variance(a,false));
+  double r = 0;
+  for(double v:a) r += pow( (v - mu)/sigma, 4 );
+  return r/a.size() - 3;
+}
+
+
+
+//Independent implementation using https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
+struct statsTest{
+  double n;
+  double mu;
+  double M2;
+  double M3;
+  double M4;
+
+  statsTest(): n(0), mu(0), M2(0), M3(0), M4(0){}
+  statsTest(const std::vector<double> &v){
+    n = v.size();
+    mu = 0;
+    M2 = 0;
+    M3 = 0;
+    M4 = 0;
+    for(double e: v)
+      mu += e;
+    mu /= n;
+    
+    for(double e: v){
+      M2 += pow(e - mu,2);
+      M3 += pow(e - mu,3);
+      M4 += pow(e - mu,4);
+    }
+      
+  }
+
+  double variance() const{ //includes Bessel's correction
+    return M2/(n-1.);
+  }
+  double mean() const{
+    return mu;
+  }
+  double skewness() const{
+    return M3/n/pow(M2/n,3./2.);
+  }
+  double kurtosis() const{
+    return M4/n/pow(M2/n,2) - 3;
+  }
+
+};
+
+statsTest operator+(const statsTest &a, const statsTest &b){ 
+  statsTest out;
+  out.n = a.n + b.n;
+  double delta = b.mu - a.mu;
+  out.mu = a.mu + delta * b.n/out.n;
+  out.M2 = a.M2 + b.M2 + pow(delta,2) * a.n * b.n / out.n;
+  out.M3 = a.M3 + b.M3 + pow(delta,3) * a.n * b.n * (a.n - b.n) / pow(out.n,2) + 3*delta*(a.n*b.M2 - b.n*a.M2)/out.n;
+  out.M4 = a.M4 + b.M4 + pow(delta,4) * a.n * b.n * (a.n*a.n - a.n*b.n + b.n*b.n)/pow(out.n,3) + 6*pow(delta,2) * (a.n*a.n*b.M2 + b.n*b.n*a.M2 ) / pow(out.n,2)
+    + 4*delta*( a.n*b.M3 - b.n*a.M3 ) / out.n;
+  
+  return out;
+}
+
+bool compare(const statsTest &a, const statsTest &b, const double tol = 1e-12){
+  bool ret = true;
+#define COM(A) if(2.*fabs( a. A - b. A )/(a. A  + b. A) > tol){ std::cout << #A << " " << a. A << " " << b. A << std::endl; ret = false; }
+  COM(n);
+  COM(mu);
+  COM(M2);
+  COM(M3);
+  COM(M4);
+  return ret;
+#undef COM
+}
+
+
+bool compare(const statsTest &a, const RunStats &b, const double tol = 1e-12){
+  const RunStats::State &sb = b.get_state();
+
+  // double eta; /**< mean */
+  // double rho; /**< = M2 = \sum_i (x_i - \bar x)^2 */
+  // double tau; /**< = M3 = \sum_i (x_i - \bar x)^3 */
+  // double phi; /**< = M4 = \sum_i (x_i - \bar x)^4 */
+
+  bool ret = true;
+#define COM(A,B) if( \
+		    (fabs(a. A)<=tol && fabs(sb. B)>tol) || \
+		    (fabs(sb. B)<=tol && fabs(a. A)>tol) ||		\
+		    (fabs(a. A)>tol && fabs(sb. B) > tol && 2.*fabs( a. A - sb. B )/(a. A  + sb. B) > tol) \
+		     ){ std::cout << #A << " " << a. A << " " << #B << " " << sb. B << std::endl; ret = false; }
+  COM(n,count);
+  COM(mu,eta);
+  COM(M2,rho);
+  COM(M3,tau);
+  COM(M4,phi);
+  return ret;
+#undef COM
+}
+
+
+TEST(TestRunStats, TestIndependentImplementation){
+  //Test that summing two RunStats is the same as if the data were collected by a single RunStats instance
+  std::vector<std::vector<double> > all_vals = {
+    {160,150,140,122,103,77,33,22,19,7,1},
+    {77,33,22,19,7,1},
+    {77,33,22,19},
+    {-0.2, -0.5, 0.7, -0.4},
+    {3.14,6.28,9.99,10.22},
+    {1000,2000,3000,4000},
+    {22,-22,22,-22}
+  };
+  for(auto const &vals: all_vals){
+    std::vector<double> data_a, data_b;
+
+    int na = vals.size() / 2;
+    int nb = vals.size() - na;
+    for(int i=0;i<na;i++)
+      data_a.push_back(vals[i]);
+    for(int i=na;i<na+nb;i++)
+      data_b.push_back(vals[i]);
+  
+    for(int i=0;i<vals.size();i++){
+      std::cout << vals[i] << " ";
+    }
+    std::cout << std::endl;
+
+    statsTest a(data_a), b(data_b), c(vals);
+
+    ASSERT_NEAR(c.mean(), mean(vals), 1e-3);
+    ASSERT_NEAR(c.variance(), variance(vals), 1e-3);
+    ASSERT_NEAR(c.skewness(), skewness(vals), 1e-3);
+    ASSERT_NEAR(c.kurtosis(), kurtosis(vals), 1e-3);
+    std::cout << "Full dist mean " << c.mean() << " var " << c.variance() << " skewness " << c.skewness() << " kurtosis " << c.kurtosis() << " match expected" << std::endl;
+    
+    statsTest sum = a + b;
+
+    bool result = compare(c, sum, 1e-10);
+
+    std::cout << "Result a+b: " << (result?"pass":"fail") << std::endl;
+
+    EXPECT_EQ(result,true);
+    ASSERT_NEAR(c.mean(), sum.mean(),1e-5);
+    ASSERT_NEAR(c.variance(), sum.variance(),1e-5);
+    ASSERT_NEAR(c.skewness(), sum.skewness(),1e-5);
+    ASSERT_NEAR(c.kurtosis(), sum.kurtosis(),1e-5);
+  }
+
+}
+ 
+
+TEST(TestRunStats, TestSumCombine){
+  //Test that summing two RunStats is the same as if the data were collected by a single RunStats instance
+  std::vector<std::vector<double> > all_vals = {
+    {160,150,140,122,103,77,33,22,19,7,1},
+    {77,33,22,19,7,1},
+    {77,33,22,19},
+    {-0.2, -0.5, 0.7, -0.4},
+    {3.14,6.28,9.99,10.22},
+    {1000,2000,3000,4000},
+    {22,-22,22,-22}
+  };
+  for(auto const &vals: all_vals){
+    RunStats a(true),b(true),c(true);
+
+    int na = vals.size() / 2;
+    int nb = vals.size() - na;
+    for(int i=0;i<na;i++)
+      a.push(vals[i]);
+    for(int i=na;i<na+nb;i++)
+      b.push(vals[i]);
+
+    for(int i=0;i<vals.size();i++){
+      std::cout << vals[i] << " ";
+      c.push(vals[i]);
+    }
+    std::cout << std::endl;
+
+    //Check against independent implementation
+    std::vector<double> data_a, data_b;
+    for(int i=0;i<na;i++)
+      data_a.push_back(vals[i]);
+    for(int i=na;i<na+nb;i++)
+      data_b.push_back(vals[i]);
+    
+    statsTest ia(data_a), ib(data_b), ic(vals);
+    
+    std::cout << "Comparing distribution 'a' to independent implementation" << std::endl;
+    bool result = compare(ia,a, 1e-10);
+    std::cout << "Result: " << (result?"pass":"fail") << std::endl;
+    ASSERT_EQ(result,true);	
+
+    std::cout << "Comparing distribution 'b' to independent implementation" << std::endl;
+    result = compare(ib,b, 1e-10);
+    std::cout << "Result: " << (result?"pass":"fail") << std::endl;   
+    ASSERT_EQ(result,true);
+
+    std::cout << "Comparing distribution 'c' to independent implementation" << std::endl;
+    result = compare(ic,c, 1e-10);
+    std::cout << "Result: " << (result?"pass":"fail") << std::endl;
+    ASSERT_EQ(result,true);    
+
+    std::cout << "Comparing distribution 'c' moments to textbook definitions" << std::endl;
+    ASSERT_NEAR(c.mean(), mean(vals), 1e-3);
+    ASSERT_NEAR(c.variance(), variance(vals), 1e-3);
+    ASSERT_NEAR(c.skewness(), skewness(vals), 1e-3);
+    ASSERT_NEAR(c.kurtosis(), kurtosis(vals), 1e-3);
+    std::cout << "Full dist mean " << c.mean() << " var " << c.variance() << " skewness " << c.skewness() << " kurtosis " << c.kurtosis() << " match expected" << std::endl;
+
+    statsTest isum = ia+ib;
+    ASSERT_EQ(compare(ic,isum,1e-10),true);
+    
+    RunStats sum = a + b;
+    std::cout << "Comparing combined distribution 'a+b' to independent implementation" << std::endl;
+    result = compare(ic,sum, 1e-10);
+    std::cout << "Result: " << (result?"pass":"fail") << std::endl;
+    ASSERT_EQ(result,true);
+
+
+    std::cout << "Comparing combined distribution 'a+b' to 'c'" << std::endl;
+    result = compare(c, sum, 1e-10);
+    std::cout << "Result: " << (result?"pass":"fail") << std::endl;
+    EXPECT_EQ(result,true);
+
+    sum = b + a;
+    
+    std::cout << "Comparing combined distribution 'b+a' to 'c'" << std::endl;
+    result = compare(c, sum, 1e-10);
+    std::cout << "Result: " << (result?"pass":"fail") << std::endl;
+    EXPECT_EQ(result,true);
+  }
+}
+ 
+
 TEST(TestRunStats, TestStateToFromJSON){
   RunStats stats;
   for(int i=0;i<100;i++) stats.push(i);
@@ -53,7 +310,6 @@ TEST(TestRunStats, serialize){
   EXPECT_EQ(stats, stats_rd);
 }
 
-