Ignore:
Timestamp:
Sep 24, 2011, 10:53:04 PM (10 years ago)
Author:
Peter
Message:

refs #671. Follow Knuth's algorithm to calculate average and variance online

File:
1 edited

Legend:

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

    r2119 r2558  
    66  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér
    77  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
     8  Copyright (C) 2011 Peter Johansson
    89
    910  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    3334
    3435  Averager::Averager(void)
    35     : n_(0), x_(0), xx_(0)
     36    : n_(0), mean_(0), m2_(0)
    3637  {
    3738  }
    3839
     40
    3941  Averager::Averager(double x, double xx, long n)
    40     : n_(n), x_(x), xx_(xx)
     42    : n_(n), mean_(x/n), m2_(xx-x*x/n)
    4143  {
    4244  }
    4345
     46
    4447  Averager::Averager(const Averager& a)
    45     : n_(a.n_), x_(a.x_), xx_(a.xx_)
     48    : n_(a.n_), mean_(a.mean_), m2_(a.m2_)
    4649  {
    4750  }
    4851
     52
    4953  void Averager::add(double d, long n)
    5054  {
     55    double delta = d - mean_;
     56    mean_ += n*delta/(n+n_);
     57    m2_ += delta*delta*n_*n/(n_+n);
    5158    n_  += n;
    5259    assert(n_>-1);
    53     x_  += n*d;
    54     xx_ += n*d*d;
    5560  }
    5661
     
    6267  double Averager::mean(void) const
    6368  {
    64     return x_/n_;
     69    return mean_;
    6570  }
    6671
     
    7277  void Averager::rescale(double a)
    7378  {
    74     x_  *= a;
    75     xx_ *= a*a;
     79    mean_  *= a;
     80    m2_ *= a*a;
    7681  }
    7782
     
    7984  {
    8085    n_=0;
    81     x_=xx_=0.0;
     86    mean_=m2_=0.0;
    8287  }
    8388
     
    99104  double Averager::sum_x(void)  const
    100105  {
    101     return x_;
     106    return n_*mean_;
    102107  }
    103108
    104109  double Averager::sum_xx(void) const
    105110  {
    106     return xx_;
     111    return m2_+n_*mean_*mean_;
    107112  }
    108113
    109114  double Averager::sum_xx_centered(void)  const
    110115  {
    111     return xx_-x_*x_/n_;
     116    return m2_;
    112117  }
    113118
    114119  double Averager::variance(double m) const
    115120  {
    116     return (xx_ - 2*m*x_ + m*m*n()) /n_;
     121    return sum_xx()/n() + m*m - 2*m*mean();
    117122  }
    118123
     
    130135  const Averager& Averager::operator=(const Averager& a)
    131136  {
    132     n_  = a.n_;
    133     x_  = a.x_;
    134     xx_ = a.xx_;
     137    if (this != &a) { // avoid self-assignment
     138      n_  = a.n_;
     139      mean_  = a.mean_;
     140      m2_ = a.m2_;
     141    }
    135142    return *this;
    136143  }
     
    138145  const Averager& Averager::operator+=(const Averager& a)
    139146  {
     147    mean_ += (n()*mean() + a.n()*a.mean()) / (n() + a.n());
     148    double delta = mean_ = a.mean();
     149    m2_ += a.m2_ + n()*a.n()*delta*delta/(n()+a.n());
    140150    n_+=a.n_;
    141     x_+=a.x_;
    142     xx_+=a.xx_;
    143151    return *this;
    144152  }
Note: See TracChangeset for help on using the changeset viewer.