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

#ifndef _theplu_yat_statistics_averagerweighted_
#define _theplu_yat_statistics_averagerweighted_

// $Id: AveragerWeighted.h 756 2007-02-19 18:14:23Z 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 <cassert>
#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:

    ///
    /// @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);

    ///
    /// 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);

    ///
    /// @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_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;
    
    ///
    ///  @return \f$ \sum w_i x_i^2 \f$
    ///
    double sum_wxx(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_;
  };

  // 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