# Changeset 486

Ignore:
Timestamp:
Jan 4, 2006, 5:41:48 PM (16 years ago)
Message:

aligned AveragerWeighted? with Averager and added tests for checking this see class doc

Location:
trunk
Files:
4 edited

Unmodified
Removed
• ## trunk/lib/statistics/AveragerWeighted.cc

 r420 #include #include #include #include namespace statistics{ void AveragerWeighted::add(const gslapi::vector& x, const gslapi::vector& w) { assert(x.size()==w.size()); for (size_t i=0; i
• ## trunk/lib/statistics/AveragerWeighted.h

 r483 #include #include //#include namespace theplu{ class gslapi::vector; namespace statistics{ /// /// Class to calulate simple (first and second moments) averages /// with weights. /// @brief Class to calulate averages with weights. /// /// @see Averager /// There are several different reasons why a statistical analysis /// needs to adjust for weighting. In the litterature reasons are /// mainly divided into two kinds of weights - probablity weights /// and analytical weights. 1) Analytical weights are appropriate /// for scientific experiments where some measurements are known to /// be more precise than others. The larger weight a measurement has /// the more precise is is assumed to be, or more formally the /// weight is proportional to the reciprocal variance /// \f$\sigma_i^2 = \frac{\sigma^2}{w_i} \f$. 2) Probablity weights /// are used for the situation when calculating averages over a /// distribution \f$f \f$ , but sampling from a distribution \f$f' /// \f$. Compensating for this discrepancy averages of observables /// are taken to be \f$\sum \frac{f}{f'}X \f$ For further discussion: /// Weighted Statistics document
/// /// If nothing else stated, each function fulfills the /// following:
• Setting a weight to zero corresponds to /// removing the data point from the dataset.
• Setting all /// weights to unity, the yields the same result as from /// corresponding function in Averager.
• Rescaling weights /// does not change the performance of the object.
/// /// @see Averager AveragerPair AveragerPairWeighted /// class AveragerWeighted /// /// Adding each value in vector \a x and corresponding value in /// weight vector \a w. /// void add(const gslapi::vector& x, const gslapi::vector& w); /// /// Calculating the weighted mean /// /// /// The standard error is calculated as \f$\sqrt{\frac{\sum /// w_i^2}{(\sum w_i)^2-\sum w_i^2}\frac{\sum w_i(x_i-m)^2}{\sum /// w_i}} \f$ /// /// @return standard deviation of mean /// inline double standard_error(void)  const { return sqrt(sum_ww()/((sum_w()*sum_w())-sum_ww()) * sum_xx_centered()/sum_w()); } /// /// The standard deviation is defined as the square root of the /// variance. /// variance(). /// /// @return The standard deviation, root of the variance(). /// /// Calculating the sum of weights: \f$\sum /// w_i \f$ @return sum of weights /// Calculates standard deviation of the mean(). Variance from the /// weights are here neglected. This is true when the weight is /// known before the measurement. In case this is not a good /// approximation, use bootstrapping to estimate the error. /// /// @return \f$\frac{\sum w^2}{\left(\sum w\right)^3}\sum w(x-m)^2 \f$ /// where \f$m \f$ is the mean() /// inline double standard_error(void)  const { return sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered()); } /// /// Calculating the sum of weights /// /// @return \f$\sum w_i \f$ /// inline double sum_w(void) const /// /// @return \f$\sum_i w_i (x_i-m)^2\f$ /// inline double sum_xx_centered(void) const { return sum_wxx() - mean()*mean()*sum_w(); } /// /// The variance is calculated as \f$\frac{\sum w_i (x_i - m)^2 /// }{\sum w_i} \f$. /// }{\sum w_i} \f$, where \a m is the known mean. /// /// @return Variance when the mean is known to be \a m. /// inline double variance(const double m) const { return (sum_wxx()-2*m*sum_wx())/sum_w()+m*m; } /// /// The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2 /// }{\sum w_i} \f$, where \a m is the mean(). Here the weight are /// interpreted as probability weights. For analytical weights the /// variance has no meaning as each data point has its own /// variance. /// /// @return The variance. inline double variance(void) const { return sum_xx_centered()/sum_w(); } /// /// The variance is calculated as \f$ \frac{\sum w_i}{\left(\sum /// w_i\right)^2 - \sum w_i^2} \sum w_i (x_i - m)^2 \f$. This /// variance is the same as the one implemented in GSL. /// /// @return The variance. /// inline double variance_unbiased(void) const { return sum_ww() / (sum_w()*sum_w()-sum_ww()) * sum_xx_centered(); } /// /// @return \f$ \sum_i (x_i-m)^2\f$/// /// @note this function is likely to be changed, when this class /// will be revised. /// inline double sum_xsqr_centered(void) const { return sum_xx_centered(); } /// /// @return \f$ \sum_i w_i (x_i-m)^2\f\$ /// inline double sum_xx_centered(void) const { return sum_wxx() - mean()*mean()*sum_w(); } /// /// operator to add a AveragerWeighted /// /// The AveragerWeighted output operator /// std::ostream& operator<<(std::ostream& s,const AveragerWeighted&); ///std::ostream& operator<<(std::ostream& s,const AveragerWeighted&); }} // of namespace statistics and namespace theplu
• ## trunk/lib/statistics/tScore.cc

 r483 } double diff = positive.mean() - negative.mean(); double s=sqrt((positive.sum_xsqr_centered()+negative.sum_xsqr_centered())/ double s=sqrt((positive.sum_xx_centered()+negative.sum_xx_centered())/ (positive.sum_w()+negative.sum_w()-2)); t_=diff/s/(1.0/positive.sum_w()+1.0/negative.sum_w());
• ## trunk/test/averager_test.cc

 r484 #include #include #include #include bool equal(const Averager& a, const Averager& b) { return a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance(); return (a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance()); } bool equal(const AveragerWeighted& a, const AveragerWeighted& b, double tol, std::ostream* error) { bool equal = true; if  ( fabs(a.mean()-b.mean())>tol){ equal=false; *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl; } if ( fabs(a.variance()-b.variance())>tol ) { equal=false; *error << "error for variance:\t" << a.variance() << " " << b.variance() << std::endl; } if ( fabs(a.standard_error()-b.standard_error())>tol ) { equal =false; *error << "error for standard error:\t" << std::endl; } return equal; } bool ok = true; // Testing Averager *error << "testing Averager" << std::endl; Averager a; } // Testing AveragerWeighted *error << "testing AveragerWeighted" << std::endl; theplu::gslapi::vector x(3,0); x(0)=0; x(1)=1; x(2)=2; theplu::gslapi::vector w(3,1); theplu::statistics::AveragerWeighted aw; aw.add(x,w); a.reset(); a.add(x); double tol=0.0001; if (fabs(a.mean()-aw.mean())>tol || fabs(a.variance()-aw.variance())>tol || fabs(a.standard_error()-aw.standard_error())>tol || fabs(a.std()-aw.std())>tol){ *error << "error: AveragerWeighted with unitary weights should " << "be equal to Averager" << std::endl; *error << "\tAverager\tAveragerWeighted" << std::endl; *error << "mean:\t" << a.mean() << "\t" << aw.mean() << std::endl; *error << "var:\t" << a.variance() << "\t" << aw.variance() << std::endl; *error << "ste:\t" << a.standard_error() << "\t" << aw.standard_error() << std::endl; *error << "std:\t" << a.std() << "\t" << aw.std() << std::endl; ok=false; } AveragerWeighted* aw2 = new AveragerWeighted(aw); if (!equal(aw,*aw2,tol,error)){ *error << "error: AveragerWeighted copy constructor " << std::endl; ok=false; } aw2->add(12,0); if (!equal(aw,*aw2,tol,error)){ *error << "error: AveragerWeighted adding a data point with weight=0 " << "should make no change " << std::endl; ok=false; } aw2->reset(); w.scale(17); aw2->add(x,w); if (!equal(aw,*aw2,tol,error)){ *error << "error: AveragerWeighted rescaling weights " << "should make no change " << std::endl; ok=false; } delete aw2; theplu::statistics::AveragerPair ap; for (int i=0; i<10; i++)
Note: See TracChangeset for help on using the changeset viewer.