Changeset 2558


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

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/tukey.cc

    r2510 r2558  
    6666  TukeyBiweightEstimator estimator(5, false);
    6767  double x = estimator(data.begin(), data.end());
    68   suite.add(suite.equal(x, correct));
     68  suite.add(suite.equal(x, correct, 2));
    6969  double y = estimator(data.begin(), data.end());
    7070  TukeyBiweightEstimator estimator2(5, true);
  • 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  }
  • trunk/yat/statistics/Averager.h

    r2161 r2558  
    88  Copyright (C) 2005, 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
     
    179179  private:
    180180    long  n_;
    181     double  x_, xx_;
     181    double  mean_, m2_;
    182182  };
    183183 
Note: See TracChangeset for help on using the changeset viewer.