Changeset 2562
 Timestamp:
 Sep 25, 2011, 8:16:56 PM (10 years ago)
 Location:
 trunk/yat/statistics
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/yat/statistics/AveragerWeighted.cc
r2451 r2562 22 22 23 23 #include "AveragerWeighted.h" 24 #include "Averager.h" 24 25 #include <cmath> 25 26 26 27 namespace theplu { … … 29 30 30 31 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) 32 33 { 33 34 } 34 35 36 35 37 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_) 40 40 { 41 41 } 42 42 43 43 44 void AveragerWeighted::add(const double d, const double w) … … 45 46 if (!w) 46 47 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; 48 57 } 58 49 59 50 60 double AveragerWeighted::mean(void) const 51 61 { 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_; 53 65 } 66 54 67 55 68 double AveragerWeighted::n(void) const … … 58 71 } 59 72 73 60 74 void AveragerWeighted::rescale(double a) 61 75 { 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; 65 81 } 82 66 83 67 84 void AveragerWeighted::reset(void) 68 85 { 69 wx_.reset(); w_.reset(); wwx_=0; wxx_=0; 86 AveragerWeighted tmp; 87 *this = tmp; 70 88 } 89 71 90 72 91 double AveragerWeighted::std(void) const 73 92 { 74 return s qrt(variance());93 return std::sqrt(variance()); 75 94 } 95 76 96 77 97 double AveragerWeighted::standard_error(void) const 78 98 { 79 return s qrt(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()); 80 100 } 101 81 102 82 103 double AveragerWeighted::sum_w(void) const 83 104 { 84 return w_ .sum_x();105 return w_; 85 106 } 107 86 108 87 109 double AveragerWeighted::sum_ww(void) const 88 110 { 89 return w _.sum_xx();111 return ww_; 90 112 } 113 91 114 92 115 double AveragerWeighted::sum_wwx(void) const 93 116 { 94 return wwx_;117 return m2_ + mean_*mean_*w_; 95 118 } 119 96 120 97 121 double AveragerWeighted::sum_wwxx(void) const 98 122 { 99 return w x_.sum_xx();123 return wwxx_; 100 124 } 125 101 126 102 127 double AveragerWeighted::sum_wx(void) const 103 128 { 104 return wx_.sum_x(); 129 if (!w_) 130 return 0.0; 131 return mean_ * w_; 105 132 } 133 106 134 107 135 double AveragerWeighted::sum_wxx(void) const … … 110 138 } 111 139 140 112 141 double AveragerWeighted::sum_xx_centered(void) const 113 142 { 114 return sum_wxx()  mean()*mean()*sum_w();143 return m2_; 115 144 } 145 116 146 117 147 double AveragerWeighted::variance(const double m) const … … 120 150 } 121 151 152 122 153 double AveragerWeighted::variance(void) const 123 154 { … … 125 156 } 126 157 127 const Averager& AveragerWeighted::wx(void) const 158 159 const AveragerWeighted& 160 AveragerWeighted::operator+=(const AveragerWeighted& a) 128 161 { 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_; 140 169 return *this; 141 170 } 
trunk/yat/statistics/AveragerWeighted.h
r2202 r2562 26 26 */ 27 27 28 #include "Averager.h"29 28 #include "yat/utility/iterator_traits.h" 30 29 … … 202 201 double sum_wwxx(void) const; 203 202 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_; 210 208 double wxx_; 209 210 // double wwx; // m2_ + m*m*w 211 // double wx; // w*mean_ 211 212 }; 212 213
Note: See TracChangeset
for help on using the changeset viewer.