#ifndef _theplu_yat_statistics_averagerweighted_
#define _theplu_yat_statistics_averagerweighted_
// $Id: AveragerWeighted.h 1275 2008-04-11 06:10:12Z 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://trac.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 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 "yat/utility/iterator_traits.h"
#include
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:

- 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
{
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
void add(AveragerWeighted& a, Iter first, Iter last)
{
for ( ; first != last; ++first)
a.add(utility::iterator_traits().data(first),
utility::iterator_traits().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
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