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

#ifndef _theplu_yat_statistics_averagerweighted_
#define _theplu_yat_statistics_averagerweighted_

// $Id: AveragerWeighted.h 1486 2008-09-09 21:17:19Z jari $

/*
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson, Cecilia Ritz
  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
  Copyright (C) 2008 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/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 3 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 "yat/utility/iterator_traits.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:
  /// see \ref weighted_statistics
  /// 
  /// 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:

    ///
    /// @brief The default constructor
    /// 
    AveragerWeighted(void);

    ///
    /// @brief The copy constructor
    ///
    AveragerWeighted(const AveragerWeighted&);

    ///
    /// Adding a data point \a d, with weight \a w (default is 1)
    ///
    void  add(const double d,const double w=1);

    ///
    /// @brief Calculate the weighted mean 
    ///
    /// @return \f$ \frac{\sum w_ix_i}{\sum w_i} \f$ 
    ///
    double mean(void) const;

    ///
    /// @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$
    /// 
    double n(void) const;

    ///
    /// @brief Rescale object.
    ///
    /// Each data point is rescaled as \f$ x = a * x \f$
    ///
    void rescale(double a);

    ///
    /// @brief Reset everything to zero.
    ///
    void reset(void);

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

    ///
    /// @brief 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() 
    ///
    double standard_error(void) const;

    ///
    /// Calculating the sum of weights  
    ///
    /// @return \f$ \sum w_i \f$
    ///
    double sum_w(void) const;

    ///
    /// @return \f$ \sum w_i^2 \f$
    ///
    double sum_ww(void) const;

    ///
    /// \f$ \sum w_ix_i \f$ 
    ///
    /// @return weighted sum of x
    ///
    double sum_wx(void) const;

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

    ///
    /// @return \f$ \sum_i w_i (x_i-m)^2\f$
    ///
    double sum_xx_centered(void) const;

    /**
       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.
    */
    double variance(const double m) const;

    /**
       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.
    */
    double variance(void) const;


  private:
    ///
    ///  @return \f$ \sum w_i^2x_i \f$
    ///
    double sum_wwx(void) const;

    ///
    ///  @return \f$ \sum w_i^2x_i^2 \f$
    ///
    double sum_wwxx(void) const;
    
    const Averager& wx(void) const;
    const Averager& w(void) const;

    ///
    /// operator to add a AveragerWeighted
    ///
    const AveragerWeighted& operator+=(const AveragerWeighted&);
    
    Averager w_;
    Averager wx_;
    double wwx_;
    double wxx_;
  };

  /**
     \brief adding a range of values to AveragerWeighted \a a

     If iterator is non-weighted unitary weights are used.
   */
  template <typename Iter>
  void add(AveragerWeighted& a, Iter first, Iter last)
  {
    for ( ; first != last; ++first)
      a.add(utility::iterator_traits<Iter>().data(first),
            utility::iterator_traits<Iter>().weight(first));
  }

  /**
     \brief add values from two ranges to AveragerWeighted \a a 

     Add data from range [first1, last1) with their corresponding
     weight in range [first2, first2 + distance(first, last) ).

     Requirement: Iter1 and Iter2 are unweighted iterators.
   */
  template <typename Iter1, typename Iter2>
  void add(AveragerWeighted& a, Iter1 first1, Iter1 last1, Iter2 first2)
  {
    utility::check_iterator_is_unweighted(first1);
    utility::check_iterator_is_unweighted(first2);
    for ( ; first1 != last1; ++first1, ++first2)
      a.add(*first1, *first2);
  }

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

#endif