Changeset 486
 Timestamp:
 Jan 4, 2006, 5:41:48 PM (17 years ago)
 Location:
 trunk
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

trunk/lib/statistics/AveragerWeighted.cc
r420 r486 3 3 #include <c++_tools/statistics/AveragerWeighted.h> 4 4 #include <c++_tools/statistics/Averager.h> 5 #include <c++_tools/gslapi/vector.h> 5 6 6 7 #include <ostream> … … 10 11 namespace statistics{ 11 12 13 14 void AveragerWeighted::add(const gslapi::vector& x, const gslapi::vector& w) 15 { 16 assert(x.size()==w.size()); 17 for (size_t i=0; i<x.size(); i++) 18 add(x(i),w(i)); 19 } 20 21 12 22 13 23 AveragerWeighted AveragerWeighted::operator+=(const AveragerWeighted& a) 
trunk/lib/statistics/AveragerWeighted.h
r483 r486 7 7 8 8 #include <cmath> 9 #include <ostream>9 //#include <ostream> 10 10 11 11 namespace theplu{ 12 class gslapi::vector; 13 12 14 namespace statistics{ 13 15 14 16 /// 15 /// Class to calulate simple (first and second moments) averages 16 /// with weights. 17 /// @brief Class to calulate averages with weights. 17 18 /// 18 /// @see Averager 19 /// There are several different reasons why a statistical analysis 20 /// needs to adjust for weighting. In the litterature reasons are 21 /// mainly divided into two kinds of weights  probablity weights 22 /// and analytical weights. 1) Analytical weights are appropriate 23 /// for scientific experiments where some measurements are known to 24 /// be more precise than others. The larger weight a measurement has 25 /// the more precise is is assumed to be, or more formally the 26 /// weight is proportional to the reciprocal variance 27 /// \f$ \sigma_i^2 = \frac{\sigma^2}{w_i} \f$. 2) Probablity weights 28 /// are used for the situation when calculating averages over a 29 /// distribution \f$ f \f$ , but sampling from a distribution \f$ f' 30 /// \f$. Compensating for this discrepancy averages of observables 31 /// are taken to be \f$ \sum \frac{f}{f'}X \f$ For further discussion: 32 /// <a href="../Statistics/index.html">Weighted Statistics document</a><br> 33 /// 34 /// If nothing else stated, each function fulfills the 35 /// following:<br> <ul><li>Setting a weight to zero corresponds to 36 /// removing the data point from the dataset.</li><li> Setting all 37 /// weights to unity, the yields the same result as from 38 /// corresponding function in Averager.</li><li> Rescaling weights 39 /// does not change the performance of the object.</li></ul> 40 /// 41 /// @see Averager AveragerPair AveragerPairWeighted 19 42 /// 20 43 class AveragerWeighted … … 43 66 44 67 /// 68 /// Adding each value in vector \a x and corresponding value in 69 /// weight vector \a w. 70 /// 71 void add(const gslapi::vector& x, const gslapi::vector& w); 72 73 /// 45 74 /// Calculating the weighted mean 46 75 /// … … 62 91 63 92 /// 64 /// The standard error is calculated as \f$ \sqrt{\frac{\sum65 /// w_i^2}{(\sum w_i)^2\sum w_i^2}\frac{\sum w_i(x_im)^2}{\sum66 /// w_i}} \f$67 ///68 /// @return standard deviation of mean69 ///70 inline double standard_error(void) const71 { return sqrt(sum_ww()/((sum_w()*sum_w())sum_ww()) *72 sum_xx_centered()/sum_w()); }73 74 ///75 93 /// The standard deviation is defined as the square root of the 76 /// variance .94 /// variance(). 77 95 /// 78 96 /// @return The standard deviation, root of the variance(). … … 81 99 82 100 /// 83 /// Calculating the sum of weights: \f$ \sum 84 /// w_i \f$ @return sum of weights 101 /// Calculates standard deviation of the mean(). Variance from the 102 /// weights are here neglected. This is true when the weight is 103 /// known before the measurement. In case this is not a good 104 /// approximation, use bootstrapping to estimate the error. 105 /// 106 /// @return \f$ \frac{\sum w^2}{\left(\sum w\right)^3}\sum w(xm)^2 \f$ 107 /// where \f$ m \f$ is the mean() 108 /// 109 inline double standard_error(void) const 110 { return sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * 111 sum_xx_centered()); } 112 113 /// 114 /// Calculating the sum of weights 115 /// 116 /// @return \f$ \sum w_i \f$ 85 117 /// 86 118 inline double sum_w(void) const … … 88 120 89 121 /// 122 /// @return \f$ \sum_i w_i (x_im)^2\f$ 123 /// 124 inline double sum_xx_centered(void) const 125 { return sum_wxx()  mean()*mean()*sum_w(); } 126 127 /// 90 128 /// The variance is calculated as \f$ \frac{\sum w_i (x_i  m)^2 91 /// }{\sum w_i} \f$. 129 /// }{\sum w_i} \f$, where \a m is the known mean. 130 /// 131 /// @return Variance when the mean is known to be \a m. 132 /// 133 inline double variance(const double m) const 134 { return (sum_wxx()2*m*sum_wx())/sum_w()+m*m; } 135 136 /// 137 /// The variance is calculated as \f$ \frac{\sum w_i (x_i  m)^2 138 /// }{\sum w_i} \f$, where \a m is the mean(). Here the weight are 139 /// interpreted as probability weights. For analytical weights the 140 /// variance has no meaning as each data point has its own 141 /// variance. 92 142 /// 93 143 /// @return The variance. … … 95 145 inline double variance(void) const 96 146 { return sum_xx_centered()/sum_w(); } 97 98 ///99 /// The variance is calculated as \f$ \frac{\sum w_i}{\left(\sum100 /// w_i\right)^2  \sum w_i^2} \sum w_i (x_i  m)^2 \f$. This101 /// variance is the same as the one implemented in GSL.102 ///103 /// @return The variance.104 ///105 inline double variance_unbiased(void) const106 { return sum_ww() / (sum_w()*sum_w()sum_ww()) * sum_xx_centered(); }107 108 ///109 /// @return \f$ \sum_i (x_im)^2\f$110 ///111 /// @note this function is likely to be changed, when this class112 /// will be revised.113 ///114 inline double sum_xsqr_centered(void) const { return sum_xx_centered(); }115 147 116 148 … … 143 175 144 176 /// 145 /// @return \f$ \sum_i w_i (x_im)^2\f$146 ///147 inline double sum_xx_centered(void) const148 { return sum_wxx()  mean()*mean()*sum_w(); }149 150 ///151 177 /// operator to add a AveragerWeighted 152 178 /// … … 167 193 /// The AveragerWeighted output operator 168 194 /// 169 std::ostream& operator<<(std::ostream& s,const AveragerWeighted&);195 ///std::ostream& operator<<(std::ostream& s,const AveragerWeighted&); 170 196 171 197 }} // of namespace statistics and namespace theplu 
trunk/lib/statistics/tScore.cc
r483 r486 57 57 } 58 58 double diff = positive.mean()  negative.mean(); 59 double s=sqrt((positive.sum_x sqr_centered()+negative.sum_xsqr_centered())/59 double s=sqrt((positive.sum_xx_centered()+negative.sum_xx_centered())/ 60 60 (positive.sum_w()+negative.sum_w()2)); 61 61 t_=diff/s/(1.0/positive.sum_w()+1.0/negative.sum_w()); 
trunk/test/averager_test.cc
r484 r486 4 4 #include <c++_tools/statistics/AveragerPair.h> 5 5 #include <c++_tools/statistics/AveragerWeighted.h> 6 #include <c++_tools/gslapi/vector.h> 6 7 7 8 #include <fstream> … … 12 13 bool equal(const Averager& a, const Averager& b) 13 14 { 14 return a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance(); 15 return (a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance()); 16 } 17 18 bool equal(const AveragerWeighted& a, const AveragerWeighted& b, double tol, 19 std::ostream* error) 20 { 21 bool equal = true; 22 if ( fabs(a.mean()b.mean())>tol){ 23 equal=false; 24 *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl; 25 } 26 if ( fabs(a.variance()b.variance())>tol ) { 27 equal=false; 28 *error << "error for variance:\t" << a.variance() << " " << b.variance() 29 << std::endl; 30 } 31 if ( fabs(a.standard_error()b.standard_error())>tol ) { 32 equal =false; 33 *error << "error for standard error:\t" << std::endl; 34 } 35 return equal; 15 36 } 16 37 … … 28 49 bool ok = true; 29 50 51 // Testing Averager 30 52 *error << "testing Averager" << std::endl; 31 53 Averager a; … … 74 96 } 75 97 98 // Testing AveragerWeighted 99 *error << "testing AveragerWeighted" << std::endl; 100 theplu::gslapi::vector x(3,0); 101 x(0)=0; 102 x(1)=1; 103 x(2)=2; 104 theplu::gslapi::vector w(3,1); 76 105 theplu::statistics::AveragerWeighted aw; 106 aw.add(x,w); 107 a.reset(); 108 a.add(x); 109 double tol=0.0001; 110 if (fabs(a.mean()aw.mean())>tol  111 fabs(a.variance()aw.variance())>tol  112 fabs(a.standard_error()aw.standard_error())>tol  113 fabs(a.std()aw.std())>tol){ 114 *error << "error: AveragerWeighted with unitary weights should " 115 << "be equal to Averager" << std::endl; 116 *error << "\tAverager\tAveragerWeighted" << std::endl; 117 *error << "mean:\t" << a.mean() << "\t" << aw.mean() << std::endl; 118 *error << "var:\t" << a.variance() << "\t" << aw.variance() << std::endl; 119 *error << "ste:\t" << a.standard_error() << "\t" << aw.standard_error() 120 << std::endl; 121 *error << "std:\t" << a.std() << "\t" << aw.std() << std::endl; 122 123 ok=false; 124 } 125 126 AveragerWeighted* aw2 = new AveragerWeighted(aw); 127 if (!equal(aw,*aw2,tol,error)){ 128 *error << "error: AveragerWeighted copy constructor " << std::endl; 129 ok=false; 130 } 131 132 aw2>add(12,0); 133 if (!equal(aw,*aw2,tol,error)){ 134 *error << "error: AveragerWeighted adding a data point with weight=0 " 135 << "should make no change " << std::endl; 136 ok=false; 137 } 138 139 aw2>reset(); 140 w.scale(17); 141 aw2>add(x,w); 142 if (!equal(aw,*aw2,tol,error)){ 143 *error << "error: AveragerWeighted rescaling weights " 144 << "should make no change " << std::endl; 145 ok=false; 146 } 147 delete aw2; 148 149 77 150 theplu::statistics::AveragerPair ap; 78 151 for (int i=0; i<10; i++)
Note: See TracChangeset
for help on using the changeset viewer.