Changeset 2565


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

refs #671

Location:
trunk/yat/statistics
Files:
2 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  }
  • trunk/yat/statistics/AveragerPair.h

    r2506 r2565  
    88  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
    99  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
    10   Copyright (C) 2009, 2010 Peter Johansson
     10  Copyright (C) 2009, 2010, 2011 Peter Johansson
    1111
    1212  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    149149    Averager x_;
    150150    Averager y_;
    151     double  xy_;
     151    double  xy_centered_;
    152152
     153    void xy_add(double mx, double my, double xy_centered, long n);
    153154  };
    154155
Note: See TracChangeset for help on using the changeset viewer.