Changeset 2562


Ignore:
Timestamp:
Sep 25, 2011, 8:16:56 PM (10 years ago)
Author:
Peter
Message:

refs #671. Reimplement AveragerWeighted?

Location:
trunk/yat/statistics
Files:
2 edited

Legend:

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

    r2451 r2562  
    2222
    2323#include "AveragerWeighted.h"
    24 #include "Averager.h"
     24
     25#include <cmath>
    2526
    2627namespace theplu {
     
    2930
    3031  AveragerWeighted::AveragerWeighted(void)
    31     : w_(Averager()), wx_(Averager()), wwx_(0), wxx_(0)
     32    : mean_(0), m2_(0), w_(0), ww_(0), wwxx_(0), wxx_(0)
    3233  {
    3334  }
    3435
     36
    3537  AveragerWeighted::AveragerWeighted(const AveragerWeighted& a)
    36     : w_(a.w_),
    37       wx_(a.wx_),
    38       wwx_(a.sum_wwx()),
    39       wxx_(a.sum_wxx())
     38    : mean_(a.mean_), m2_(a.m2_), w_(a.w_),
     39      ww_(a.ww_), wwxx_(a.wwxx_), wxx_(a.wxx_)
    4040  {
    4141  }
     42
    4243
    4344  void AveragerWeighted::add(const double d, const double w)
     
    4546    if (!w)
    4647      return;
    47     w_.add(w); wx_.add(w*d); wwx_+=w*w*d; wxx_+=w*d*d;
     48    double new_w = w_ + w;
     49    double delta = d - mean_;
     50    double R = delta * w / (new_w);
     51    mean_ += R;
     52    m2_ += w_ * delta * R;
     53    w_ = new_w;
     54    ww_ += w*w;
     55    wwxx_ += w*w*d*d;
     56    wxx_ += w*d*d;
    4857  }
     58
    4959
    5060  double AveragerWeighted::mean(void) const
    5161  {
    52     return sum_wx()/sum_w();
     62    if (!w_) // no data, mean is undefined
     63      return std::numeric_limits<double>::quiet_NaN();
     64    return mean_;
    5365  }
     66
    5467
    5568  double AveragerWeighted::n(void) const
     
    5871  }
    5972
     73
    6074  void AveragerWeighted::rescale(double a)
    6175  {
    62     wx_.rescale(a);
    63     wwx_*=a;
    64     wxx_*=a*a;
     76    mean_ *= a;
     77    double a2 = a*a;
     78    m2_ *= a2;
     79    wwxx_ *= a2;
     80    wxx_*=a2;
    6581  }
     82
    6683
    6784  void AveragerWeighted::reset(void)
    6885  {
    69     wx_.reset(); w_.reset(); wwx_=0; wxx_=0;
     86    AveragerWeighted tmp;
     87    *this = tmp;
    7088  }
     89
    7190
    7291  double AveragerWeighted::std(void) const
    7392  {
    74     return sqrt(variance());
     93    return std::sqrt(variance());
    7594  }
     95
    7696
    7797  double AveragerWeighted::standard_error(void) const
    7898  {
    79     return sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered());
     99    return std::sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered());
    80100  }
     101
    81102
    82103  double AveragerWeighted::sum_w(void)  const
    83104  {
    84     return w_.sum_x();
     105    return w_;
    85106  }
     107
    86108
    87109  double AveragerWeighted::sum_ww(void) const
    88110  {
    89     return w_.sum_xx();
     111    return ww_;
    90112  }
     113
    91114
    92115  double AveragerWeighted::sum_wwx(void) const
    93116  {
    94     return wwx_;
     117    return m2_ + mean_*mean_*w_;
    95118  }
     119
    96120
    97121  double AveragerWeighted::sum_wwxx(void) const
    98122  {
    99     return wx_.sum_xx();
     123    return wwxx_;
    100124  }
     125
    101126
    102127  double AveragerWeighted::sum_wx(void) const
    103128  {
    104     return wx_.sum_x();
     129    if (!w_)
     130      return 0.0;
     131    return mean_ * w_;
    105132  }
     133
    106134
    107135  double AveragerWeighted::sum_wxx(void) const
     
    110138  }
    111139
     140
    112141  double AveragerWeighted::sum_xx_centered(void) const
    113142  {
    114     return sum_wxx() - mean()*mean()*sum_w();
     143    return m2_;
    115144  }
     145
    116146
    117147  double AveragerWeighted::variance(const double m) const
     
    120150  }
    121151
     152
    122153  double AveragerWeighted::variance(void) const
    123154  {
     
    125156  }
    126157
    127   const Averager& AveragerWeighted::wx(void) const
     158
     159  const AveragerWeighted&
     160  AveragerWeighted::operator+=(const AveragerWeighted& a)
    128161  {
    129     return wx_;
    130   }
    131 
    132   const Averager& AveragerWeighted::w(void) const
    133   {
    134     return w_;
    135   }
    136 
    137   const AveragerWeighted& AveragerWeighted::operator+=(const AveragerWeighted& a)
    138   {
    139     wx_+=a.wx(); w_+=a.w(); wwx_+=a.sum_wwx(); wxx_+=a.sum_wxx();
     162    mean_ += (sum_wx() + a.sum_wx()) / (sum_w() + a.sum_w());
     163    double delta = mean() - a.mean();
     164    m2_ += a.m2_ + sum_w()*a.sum_w()*delta*delta/(sum_w()+a.sum_w());
     165    w_ += a.w_;
     166    ww_ += a.ww_;
     167    wwxx_ += a.wwxx_;
     168    wxx_ += a.wxx_;
    140169    return *this;
    141170  }
  • trunk/yat/statistics/AveragerWeighted.h

    r2202 r2562  
    2626*/
    2727
    28 #include "Averager.h"
    2928#include "yat/utility/iterator_traits.h"
    3029
     
    202201    double sum_wwxx(void) const;
    203202   
    204     const Averager& wx(void) const;
    205     const Averager& w(void) const;
    206 
    207     Averager w_;
    208     Averager wx_;
    209     double wwx_;
     203    double mean_; // wx/w
     204    double m2_; // wxx - m*m*w
     205    double w_;
     206    double ww_;
     207    double wwxx_;
    210208    double wxx_;
     209
     210    // double wwx; // m2_ + m*m*w
     211    // double wx; // w*mean_
    211212  };
    212213
Note: See TracChangeset for help on using the changeset viewer.