Changeset 2568

Timestamp:

Sep 26, 2011, 3:38:30 AM (10 years ago)

Author:

Peter

Message:

refs #671

Location:

trunk/yat/statistics

Files:

• AveragerPairWeighted.cc (modified) (6 diffs)
• AveragerPairWeighted.h (modified) (1 diff)

trunk/yat/statistics/AveragerPairWeighted.cc

@@ -35 +35 @@
 
   AveragerPairWeighted::AveragerPairWeighted(void)
-    : wxy_(0)
+    : wxy_centered_(0)
   {
   }
@@ -51 +51 @@
     assert(!std::isnan(wy) && "wy is nan");
     double w=wx*wy;
+    xy_add(x, y, 0, w);
     x_.add(x,w);
     y_.add(y,w);
-    wxy_ += w*x*y;
   }
 
@@ -71 +71 @@
   double AveragerPairWeighted::msd(void) const
   {
-    return (x_.sum_wxx()+y_.sum_wxx()-2*wxy_)/sum_w();
+    return (x_.sum_wxx()+y_.sum_wxx()-2*sum_xy())/sum_w();
   }
 
@@ -83 +83 @@
   void AveragerPairWeighted::reset(void)
   {
-    x_.reset(); y_.reset(); wxy_=0;
+    x_.reset(); y_.reset(); wxy_centered_=0;
   }
 
@@ -95 +95 @@
   double AveragerPairWeighted::sum_xy(void) const
   {
-    return wxy_;
+    return sum_xy_centered() + x_.sum_wx()*y_.mean();
   }
 
@@ -101 +101 @@
   double AveragerPairWeighted::sum_xy_centered(void) const
   {
-    return sum_xy() - x_.sum_wx()*y_.mean();
+    return wxy_centered_;
+  }
+
+
+  void AveragerPairWeighted::xy_add(double mx, double my, double wxy_centered, 
+                                    double w)
+  {
+    if (!w)
+      return;
+    /*
+      For derivation of update formula refer to 
+      http://www.janinebennett.org/index_files/ParallelStatisticsAlgorithms.pdf
+    */
+    wxy_centered_ += wxy_centered;
+    if (sum_w())
+      wxy_centered_ += 
+        w*sum_w() / (w+sum_w()) * (mx-x_.mean()) * (my-y_.mean());
   }
 

trunk/yat/statistics/AveragerPairWeighted.h

@@ -134 +134 @@
     AveragerWeighted x_; // weighted averager with w = w_x*w_y
     AveragerWeighted y_; // weighted averager with w = w_x*w_y
-    double wxy_;
+    double wxy_centered_;
+
+    void xy_add(double mx, double my, double xy_centered, double w);
   };