Ignore:
Timestamp:
Sep 26, 2011, 2:53:48 AM (10 years ago)
Author:
Peter
Message:

refs #671

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/statistics/AveragerPair.cc

    r2119 r2565  
    2525#include "Averager.h"
    2626
     27#include <limits>
    2728#include <utility>
    2829
     
    3233
    3334  AveragerPair::AveragerPair(void)
    34     : x_(Averager()), y_(Averager()), xy_(0.0)
     35    : xy_centered_(0.0)
    3536  {
    3637  }
    3738
     39
    3840  AveragerPair::AveragerPair(const AveragerPair& a)
    39     : x_(a.x_averager()), y_(a.y_averager()), xy_(a.sum_xy())
     41    : x_(a.x_), y_(a.y_), xy_centered_(a.xy_centered_)
    4042  {
    4143  }
    4244
     45
    4346  void AveragerPair::add(const double x, const double y, const long n)
    4447  {
    45     x_.add(x,n); y_.add(y,n); xy_ += n*x*y;
     48    xy_add(x, y, 0, n);
     49    x_.add(x,n);
     50    y_.add(y,n);
    4651  }
     52
    4753
    4854  double AveragerPair::ccc(void) const
     
    5359  }
    5460
     61
    5562  double AveragerPair::correlation(void) const
    56   { return covariance() / std::sqrt(x_.variance()*y_.variance());
     63  {
     64    return covariance() / std::sqrt(x_.variance()*y_.variance());
    5765  }
     66
    5867
    5968  double AveragerPair::covariance(void) const
    6069  {
    61     return (xy_ - x_.sum_x()*y_.mean()) / n();
     70    return xy_centered_ / n();
    6271  }
     72
    6373
    6474  double AveragerPair::mean_xy(void) const
    6575  {
    66     return xy_/n();
     76    if (!n())
     77      return std::numeric_limits<double>::quiet_NaN();
     78    return sum_xy()/n();
    6779  }
     80
    6881
    6982  double AveragerPair::msd(void) const
     
    7285  }
    7386
     87
    7488  long AveragerPair::n(void) const
    7589  {
     
    7791  }
    7892
     93
    7994  void AveragerPair::reset(void)
    8095  {
    81     x_.reset(); y_.reset(); xy_=0.0;
     96    x_.reset();
     97    y_.reset();
     98    xy_centered_ = 0.0;
    8299  }
     100
    83101
    84102  const AveragerPair& AveragerPair::operator=(const AveragerPair& a)
    85103  {
    86     x_=a.x_; y_=a.y_; xy_=a.xy_;
     104    x_=a.x_; y_=a.y_; xy_centered_=a.xy_centered_;
    87105    return *this;
    88106  }
     107
    89108
    90109  double AveragerPair::sum_squared_deviation(void) const
     
    93112  }
    94113
     114
    95115  double AveragerPair::sum_xy(void) const
    96116  {
    97     return xy_;
     117    return xy_centered_ + x_.sum_x()*y_.mean();;
    98118  }
     119
    99120
    100121  double AveragerPair::sum_xy_centered(void) const
    101122  {
    102     return xy_-x_.sum_x()*y_.mean();
     123    return xy_centered_;
    103124  }
     125
     126
     127  void AveragerPair::xy_add(double mx, double my, double xy_centered, long n)
     128  {
     129    if (!n)
     130      return;
     131    /*
     132      For derivation of update formula refer to
     133      http://www.janinebennett.org/index_files/ParallelStatisticsAlgorithms.pdf
     134    */
     135    xy_centered_ += xy_centered;
     136    if (this->n())
     137      xy_centered_ += static_cast<double>(n*this->n()) / (n+this->n()) *
     138        (mx-x_.mean()) * (my-y_.mean());
     139
     140  }
     141
    104142
    105143  const Averager& AveragerPair::x_averager(void) const
     
    108146  }
    109147
     148
    110149  const Averager& AveragerPair::y_averager(void) const
    111150  {
     
    113152  }
    114153
     154
    115155  const AveragerPair& AveragerPair::operator+=(const AveragerPair& a)
    116156  {
     157    xy_add(a.x_averager().mean(), a.y_averager().mean(), a.xy_centered_,a.n());
    117158    x_+=a.x_averager();
    118159    y_+=a.y_averager();
    119     xy_+=a.sum_xy();
    120160    return *this;
    121161  }
Note: See TracChangeset for help on using the changeset viewer.