source: trunk/yat/statistics/AveragerWeighted.h @ 683

#ifndef _theplu_yat_statistics_averagerweighted_
#define _theplu_yat_statistics_averagerweighted_

// $Id: AveragerWeighted.h 683 2006-10-11 22:20:36Z jari $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "Averager.h"

#include <cmath>

namespace theplu{
namespace yat{
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) Probability 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 an array \a x and corresponding value in
    /// weight array \a w.
    ///
    /// The requirements for the types T1 and T2 of the arrays \a x
    /// and \a w are: operator[] returning an element and function 
    /// size() returning the number of elements.
    ///
    template <typename T1, typename T2>
    void add_values(const T1& x, const T2& 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; }
  
    ///
    /// @brief Weighted version of number of data points. If all
    /// weights are equal, the unweighted version is identical to the
    /// non-weighted version. Adding a data point with zero weight
    /// does not change n(). The calculated value is always smaller
    /// than the actual number of data points added to object.
    ///
    /// @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_;}


  };

  // Template implementations
  template <typename T1, typename T2>
  void AveragerWeighted::add_values(const T1& x, const T2& w)
  {
    assert(x.size()==w.size());
    for (size_t i=0; i<x.size(); i++) 
      add(x[i],w[i]);
  }

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

}}} // of namespace statistics, yat, and theplu

#endif