Changeset 2558

Timestamp:

Sep 24, 2011, 10:53:04 PM (11 years ago)

Author:

Peter

Message:

refs #671. Follow Knuth's algorithm to calculate average and variance online

Location:

trunk

Files:

• test/tukey.cc (modified) (1 diff)
• yat/statistics/Averager.cc (modified) (8 diffs)
• yat/statistics/Averager.h (modified) (2 diffs)

trunk/test/tukey.cc

@@ -66 +66 @@
   TukeyBiweightEstimator estimator(5, false);
   double x = estimator(data.begin(), data.end());
-  suite.add(suite.equal(x, correct));
+  suite.add(suite.equal(x, correct, 2));
   double y = estimator(data.begin(), data.end());
   TukeyBiweightEstimator estimator2(5, true);

trunk/yat/statistics/Averager.cc

@@ -6 +6 @@
   Copyright (C) 2006 Jari Häkkinen, Markus Ringnér
   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
+  Copyright (C) 2011 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -33 +34 @@
 
   Averager::Averager(void)
-    : n_(0), x_(0), xx_(0) 
+    : n_(0), mean_(0), m2_(0) 
   {
   }
 
+
   Averager::Averager(double x, double xx, long n)
-    : n_(n), x_(x), xx_(xx) 
+    : n_(n), mean_(x/n), m2_(xx-x*x/n) 
   {
   }
 
+
   Averager::Averager(const Averager& a)
-    : n_(a.n_), x_(a.x_), xx_(a.xx_)
+    : n_(a.n_), mean_(a.mean_), m2_(a.m2_)
   {
   }
 
+
   void Averager::add(double d, long n)
   {
+    double delta = d - mean_;
+    mean_ += n*delta/(n+n_);
+    m2_ += delta*delta*n_*n/(n_+n);
     n_  += n;
     assert(n_>-1);
-    x_  += n*d;
-    xx_ += n*d*d;
   }
 
@@ -62 +67 @@
   double Averager::mean(void) const
   {
-    return x_/n_;
+    return mean_;
   }
 
@@ -72 +77 @@
   void Averager::rescale(double a)
   {
-    x_  *= a;
-    xx_ *= a*a;
+    mean_  *= a;
+    m2_ *= a*a;
   }
 
@@ -79 +84 @@
   {
     n_=0;
-    x_=xx_=0.0;
+    mean_=m2_=0.0;
   }
 
@@ -99 +104 @@
   double Averager::sum_x(void)  const
   {
-    return x_;
+    return n_*mean_;
   }
 
   double Averager::sum_xx(void) const
   {
-    return xx_;
+    return m2_+n_*mean_*mean_;
   }
 
   double Averager::sum_xx_centered(void)  const
   {
-    return xx_-x_*x_/n_;
+    return m2_;
   }
 
   double Averager::variance(double m) const
   {
-    return (xx_ - 2*m*x_ + m*m*n()) /n_;
+    return sum_xx()/n() + m*m - 2*m*mean();
   }
 
@@ -130 +135 @@
   const Averager& Averager::operator=(const Averager& a)
   {
-    n_  = a.n_;
-    x_  = a.x_;
-    xx_ = a.xx_;
+    if (this != &a) { // avoid self-assignment
+      n_  = a.n_;
+      mean_  = a.mean_;
+      m2_ = a.m2_;
+    }
     return *this;
   }
@@ -138 +145 @@
   const Averager& Averager::operator+=(const Averager& a)
   {
+    mean_ += (n()*mean() + a.n()*a.mean()) / (n() + a.n());
+    double delta = mean_ = a.mean();
+    m2_ += a.m2_ + n()*a.n()*delta*delta/(n()+a.n());
     n_+=a.n_;
-    x_+=a.x_;
-    xx_+=a.xx_;
     return *this;
   }

trunk/yat/statistics/Averager.h

@@ -8 +8 @@
   Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
-  Copyright (C) 2009, 2010 Peter Johansson
+  Copyright (C) 2009, 2010, 2011 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -179 +179 @@
   private:
     long  n_;
-    double  x_, xx_;
+    double  mean_, m2_;
   };