Changeset 2562

Timestamp:

Sep 25, 2011, 8:16:56 PM (10 years ago)

Author:

Peter

Message:

refs #671. Reimplement AveragerWeighted?

Location:

trunk/yat/statistics

Files:

• AveragerWeighted.cc (modified) (7 diffs)
• AveragerWeighted.h (modified) (2 diffs)

trunk/yat/statistics/AveragerWeighted.cc

@@ -22 +22 @@
 
 #include "AveragerWeighted.h"
-#include "Averager.h"
+
+#include <cmath>
 
 namespace theplu {
@@ -29 +30 @@
 
   AveragerWeighted::AveragerWeighted(void)
-    : w_(Averager()), wx_(Averager()), wwx_(0), wxx_(0)
+    : mean_(0), m2_(0), w_(0), ww_(0), wwxx_(0), wxx_(0)
   {
   }
 
+
   AveragerWeighted::AveragerWeighted(const AveragerWeighted& a)
-    : w_(a.w_),
-      wx_(a.wx_),
-      wwx_(a.sum_wwx()),
-      wxx_(a.sum_wxx())
+    : mean_(a.mean_), m2_(a.m2_), w_(a.w_),
+      ww_(a.ww_), wwxx_(a.wwxx_), wxx_(a.wxx_)
   {
   }
+
 
   void AveragerWeighted::add(const double d, const double w)
@@ -45 +46 @@
     if (!w)
       return;
-    w_.add(w); wx_.add(w*d); wwx_+=w*w*d; wxx_+=w*d*d;
+    double new_w = w_ + w;
+    double delta = d - mean_;
+    double R = delta * w / (new_w);
+    mean_ += R;
+    m2_ += w_ * delta * R;
+    w_ = new_w;
+    ww_ += w*w;
+    wwxx_ += w*w*d*d; 
+    wxx_ += w*d*d;
   }
+
 
   double AveragerWeighted::mean(void) const
   {
-    return sum_wx()/sum_w();
+    if (!w_) // no data, mean is undefined
+      return std::numeric_limits<double>::quiet_NaN();
+    return mean_;
   }
+
 
   double AveragerWeighted::n(void) const
@@ -58 +71 @@
   }
 
+
   void AveragerWeighted::rescale(double a)
   {
-    wx_.rescale(a);
-    wwx_*=a;
-    wxx_*=a*a;
+    mean_ *= a;
+    double a2 = a*a;
+    m2_ *= a2;
+    wwxx_ *= a2;
+    wxx_*=a2;
   }
+
 
   void AveragerWeighted::reset(void)
   {
-    wx_.reset(); w_.reset(); wwx_=0; wxx_=0;
+    AveragerWeighted tmp;
+    *this = tmp;
   }
+
 
   double AveragerWeighted::std(void) const
   {
-    return sqrt(variance());
+    return std::sqrt(variance());
   }
+
 
   double AveragerWeighted::standard_error(void) const
   {
-    return sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered());
+    return std::sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered());
   }
+
 
   double AveragerWeighted::sum_w(void)  const
   {
-    return w_.sum_x();
+    return w_;
   }
+
 
   double AveragerWeighted::sum_ww(void) const
   {
-    return w_.sum_xx();
+    return ww_;
   }
+
 
   double AveragerWeighted::sum_wwx(void) const
   {
-    return wwx_;
+    return m2_ + mean_*mean_*w_;
   }
+
 
   double AveragerWeighted::sum_wwxx(void) const
   {
-    return wx_.sum_xx();
+    return wwxx_;
   }
+
 
   double AveragerWeighted::sum_wx(void) const
   {
-    return wx_.sum_x();
+    if (!w_)
+      return 0.0;
+    return mean_ * w_;
   }
+
 
   double AveragerWeighted::sum_wxx(void) const
@@ -110 +138 @@
   }
 
+
   double AveragerWeighted::sum_xx_centered(void) const
   {
-    return sum_wxx() - mean()*mean()*sum_w();
+    return m2_;
   }
+
 
   double AveragerWeighted::variance(const double m) const
@@ -120 +150 @@
   }
 
+
   double AveragerWeighted::variance(void) const
   {
@@ -125 +156 @@
   }
 
-  const Averager& AveragerWeighted::wx(void) const
+
+  const AveragerWeighted& 
+  AveragerWeighted::operator+=(const AveragerWeighted& a)
   {
-    return wx_;
-  }
-
-  const Averager& AveragerWeighted::w(void) const
-  {
-    return w_;
-  }
-
-  const AveragerWeighted& AveragerWeighted::operator+=(const AveragerWeighted& a)
-  {
-    wx_+=a.wx(); w_+=a.w(); wwx_+=a.sum_wwx(); wxx_+=a.sum_wxx();
+    mean_ += (sum_wx() + a.sum_wx()) / (sum_w() + a.sum_w());
+    double delta = mean() - a.mean();
+    m2_ += a.m2_ + sum_w()*a.sum_w()*delta*delta/(sum_w()+a.sum_w());
+    w_ += a.w_;
+    ww_ += a.ww_;
+    wwxx_ += a.wwxx_;
+    wxx_ += a.wxx_;
     return *this;
   }

trunk/yat/statistics/AveragerWeighted.h

@@ -26 +26 @@
 */
 
-#include "Averager.h"
 #include "yat/utility/iterator_traits.h"
 
@@ -202 +201 @@
     double sum_wwxx(void) const;
     
-    const Averager& wx(void) const;
-    const Averager& w(void) const;
-
-    Averager w_;
-    Averager wx_;
-    double wwx_;
+    double mean_; // wx/w
+    double m2_; // wxx - m*m*w
+    double w_;
+    double ww_;
+    double wwxx_;
     double wxx_;
+
+    // double wwx; // m2_ + m*m*w
+    // double wx; // w*mean_
   };