Changeset 2565

Timestamp:

Sep 26, 2011, 2:53:48 AM (11 years ago)

Author:

Peter

Message:

refs #671

Location:

trunk/yat/statistics

Files:

• AveragerPair.cc (modified) (8 diffs)
• AveragerPair.h (modified) (2 diffs)

trunk/yat/statistics/AveragerPair.cc

@@ -25 +25 @@
 #include "Averager.h"
 
+#include <limits>
 #include <utility>
 
@@ -32 +33 @@
 
   AveragerPair::AveragerPair(void)
-    : x_(Averager()), y_(Averager()), xy_(0.0)
+    : xy_centered_(0.0)
   {
   }
 
+
   AveragerPair::AveragerPair(const AveragerPair& a)
-    : x_(a.x_averager()), y_(a.y_averager()), xy_(a.sum_xy())
+    : x_(a.x_), y_(a.y_), xy_centered_(a.xy_centered_)
   {
   }
 
+
   void AveragerPair::add(const double x, const double y, const long n)
   {
-    x_.add(x,n); y_.add(y,n); xy_ += n*x*y;
+    xy_add(x, y, 0, n);
+    x_.add(x,n); 
+    y_.add(y,n); 
   }
+
 
   double AveragerPair::ccc(void) const
@@ -53 +59 @@
   }
 
+
   double AveragerPair::correlation(void) const
-  { return covariance() / std::sqrt(x_.variance()*y_.variance());
+  { 
+    return covariance() / std::sqrt(x_.variance()*y_.variance());
   }
+
 
   double AveragerPair::covariance(void) const
   {
-    return (xy_ - x_.sum_x()*y_.mean()) / n();
+    return xy_centered_ / n();
   }
+
 
   double AveragerPair::mean_xy(void) const
   {
-    return xy_/n();
+    if (!n())
+      return std::numeric_limits<double>::quiet_NaN();
+    return sum_xy()/n();
   }
+
 
   double AveragerPair::msd(void) const
@@ -72 +85 @@
   }
 
+
   long AveragerPair::n(void) const
   {
@@ -77 +91 @@
   }
 
+
   void AveragerPair::reset(void)
   {
-    x_.reset(); y_.reset(); xy_=0.0;
+    x_.reset(); 
+    y_.reset(); 
+    xy_centered_ = 0.0;
   }
+
 
   const AveragerPair& AveragerPair::operator=(const AveragerPair& a)
   {
-    x_=a.x_; y_=a.y_; xy_=a.xy_;
+    x_=a.x_; y_=a.y_; xy_centered_=a.xy_centered_;
     return *this;
   }
+
 
   double AveragerPair::sum_squared_deviation(void) const
@@ -93 +112 @@
   }
 
+
   double AveragerPair::sum_xy(void) const
   {
-    return xy_;
+    return xy_centered_ + x_.sum_x()*y_.mean();;
   }
+
 
   double AveragerPair::sum_xy_centered(void) const
   {
-    return xy_-x_.sum_x()*y_.mean();
+    return xy_centered_;
   }
+
+
+  void AveragerPair::xy_add(double mx, double my, double xy_centered, long n)
+  {
+    if (!n)
+      return;
+    /*
+      For derivation of update formula refer to 
+      http://www.janinebennett.org/index_files/ParallelStatisticsAlgorithms.pdf
+    */
+    xy_centered_ += xy_centered;
+    if (this->n())
+      xy_centered_ += static_cast<double>(n*this->n()) / (n+this->n()) *
+        (mx-x_.mean()) * (my-y_.mean());
+
+  }
+
 
   const Averager& AveragerPair::x_averager(void) const
@@ -108 +146 @@
   }
 
+
   const Averager& AveragerPair::y_averager(void) const
   {
@@ -113 +152 @@
   }
 
+
   const AveragerPair& AveragerPair::operator+=(const AveragerPair& a)
   {
+    xy_add(a.x_averager().mean(), a.y_averager().mean(), a.xy_centered_,a.n());
     x_+=a.x_averager();
     y_+=a.y_averager();
-    xy_+=a.sum_xy();
     return *this;
   }

trunk/yat/statistics/AveragerPair.h

@@ -8 +8 @@
   Copyright (C) 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
@@ -149 +149 @@
     Averager x_;
     Averager y_;
-    double  xy_;
+    double  xy_centered_;
 
+    void xy_add(double mx, double my, double xy_centered, long n);
   };