source: trunk/lib/statistics/AveragerWeighted.h @ 492

// $Id: AveragerWeighted.h 492 2006-01-09 09:38:12Z peter $

#ifndef _theplu_statistics_averager_weighted_
#define _theplu_statistics_averager_weighted_

#include <c++_tools/statistics/Averager.h>

#include <cmath>
//#include <ostream>

namespace theplu{
  class gslapi::vector;

namespace statistics{

  ///
  /// @brief Class to calulate averages with weights.
  ///
  /// 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:
  /// <a href="../Statistics/index.html">Weighted Statistics document</a><br>
  /// 
  /// If nothing else stated, each function fulfills the
  /// following:<br> <ul><li>Setting a weight to zero corresponds to
  /// removing the data point from the dataset.</li><li> Setting all
  /// weights to unity, the yields the same result as from
  /// corresponding function in Averager.</li><li> Rescaling weights
  /// does not change the performance of the object.</li></ul>
  ///
  /// @see Averager AveragerPair AveragerPairWeighted
  ///
  class AveragerWeighted
  {
  public:

    ///
    /// Default constructor
    /// 
    inline AveragerWeighted(void)
      : w_(Averager()), wx_(Averager()), wwx_(0), wxx_(0) {}

    ///
    /// Copy constructor
    ///
    inline AveragerWeighted(const AveragerWeighted& a)
      : w_(Averager(a.sum_w(),a.sum_ww(),1)),
        wx_(Averager(a.sum_wx(),a.sum_wwxx(),1)), wwx_(a.sum_wwx()),
        wxx_(a.sum_wxx()) {}

    ///
    /// adding a data point d, with weight w (default is 1)
    ///
    inline void add(const double d,const double w=1)
    { if(w==0) return; w_.add(w); wx_.add(w*d); wwx_+=w*w*d; wxx_+=w*d*d; }

    ///
    /// 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 
    ///
    /// @return \f$ \frac{\sum w_ix_i}{\sum w_i} \f$ 
    ///
    inline double mean(void) const { return sum_w() ? 
                                       sum_wx()/sum_w() : 0; }
  
    ///
    /// @return \f$ \frac{\left(\sum w_i\right)^2}{\sum w_i^2} \f$
    /// 
    inline double n(void) const { return sum_w()*sum_w()/sum_ww(); }

    ///
    /// rescale object, i.e. each data point is rescaled
    /// \f$ x = a * x \f$
    ///
    inline void rescale(double a) { wx_.rescale(a); wwx_*=a; wxx_*=a*a; }

    ///
    /// resets everything to zero
    ///
    inline void reset(void) { wx_.reset(); w_.reset(); wwx_=0; wxx_=0; }

    ///
    /// The standard deviation is defined as the square root of the
    /// variance().
    ///
    /// @return The standard deviation, root of the variance().
    ///
    inline double std(void) const { return sqrt(variance()); }

    ///
    /// 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 w_.sum_x(); }

    ///
    /// @return \f$ \sum w_i^2 \f$
    ///
    inline double sum_ww(void)  const 
    { return w_.sum_xx(); }

    ///
    /// \f$ \sum w_ix_i \f$ 
    ///
    /// @return weighted sum of x
    ///
    inline double sum_wx(void)  const 
    { return wx_.sum_x(); }

    ///
    /// @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$, 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(); }


  private:
    ///
    ///  @return \f$ \sum w_i^2x_i^2 \f$
    ///
    inline double sum_wwxx(void)  const 
    { return wx_.sum_xx(); }
    
    ///
    ///  @return \f$ \sum w_i^2x_i \f$
    ///
    inline double sum_wwx(void) const 
    { return wwx_; }

    ///
    ///  @return \f$ \sum w_i x_i^2 \f$
    ///
    inline double sum_wxx(void) const { return wxx_; }

    ///
    /// operator to add a AveragerWeighted
    ///
    AveragerWeighted operator+=(const AveragerWeighted&);
    
    Averager w_;
    Averager wx_;
    double wwx_;
    double wxx_;
    
    inline Averager wx(void) const {return wx_;}
    inline Averager w(void) const {return w_;}


  };

///
/// The AveragerWeighted output operator
///
///std::ostream& operator<<(std::ostream& s,const AveragerWeighted&);

}} // of namespace statistics and namespace theplu

#endif