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

Last change on this file since 2562 was 2562, checked in by Peter, 10 years ago

refs #671. Reimplement AveragerWeighted?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.1 KB
Line 
1#ifndef _theplu_yat_statistics_averagerweighted_
2#define _theplu_yat_statistics_averagerweighted_
3
4// $Id: AveragerWeighted.h 2562 2011-09-25 18:16:56Z peter $
5
6/*
7  Copyright (C) 2004 Jari Häkkinen, Peter Johansson, Cecilia Ritz
8  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
9  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
10  Copyright (C) 2009, 2010 Peter Johansson
11
12  This file is part of the yat library, http://dev.thep.lu.se/yat
13
14  The yat library is free software; you can redistribute it and/or
15  modify it under the terms of the GNU General Public License as
16  published by the Free Software Foundation; either version 3 of the
17  License, or (at your option) any later version.
18
19  The yat library is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22  General Public License for more details.
23
24  You should have received a copy of the GNU General Public License
25  along with yat. If not, see <http://www.gnu.org/licenses/>.
26*/
27
28#include "yat/utility/iterator_traits.h"
29
30#include <boost/concept_check.hpp>
31
32#include <cmath>
33
34namespace theplu{
35namespace yat{
36namespace statistics{
37
38  ///
39  /// @brief Class to calulate averages with weights.
40  ///
41  /// There are several different reasons why a statistical analysis
42  /// needs to adjust for weighting. In the litterature reasons are
43  /// mainly divided into two kinds of weights - probablity weights
44  /// and analytical weights. 1) Analytical weights are appropriate
45  /// for scientific experiments where some measurements are known to
46  /// be more precise than others. The larger weight a measurement has
47  /// the more precise is is assumed to be, or more formally the
48  /// weight is proportional to the reciprocal variance
49  /// \f$ \sigma_i^2 = \frac{\sigma^2}{w_i} \f$. 2) Probability weights
50  /// are used for the situation when calculating averages over a
51  /// distribution \f$ f \f$ , but sampling from a distribution \f$ f' \f$.
52  /// Compensating for this discrepancy averages of observables
53  /// are taken to be \f$ \sum \frac{f}{f'}X \f$ For further discussion:
54  /// see \ref weighted_statistics
55  ///
56  /// If nothing else stated, each function fulfills the
57  /// following:<br> <ul><li>Setting a weight to zero corresponds to
58  /// removing the data point from the dataset.</li><li> Setting all
59  /// weights to unity, the yields the same result as from
60  /// corresponding function in Averager.</li><li> Rescaling weights
61  /// does not change the performance of the object.</li></ul>
62  ///
63  /// @see Averager AveragerPair AveragerPairWeighted
64  ///
65  class AveragerWeighted
66  {
67  public:
68
69    ///
70    /// @brief The default constructor
71    ///
72    AveragerWeighted(void);
73
74    ///
75    /// @brief The copy constructor
76    ///
77    AveragerWeighted(const AveragerWeighted&);
78
79    ///
80    /// Adding a data point \a d, with weight \a w (default is 1)
81    ///
82    void  add(const double d,const double w=1);
83
84    ///
85    /// @brief Calculate the weighted mean
86    ///
87    /// @return \f$ \frac{\sum w_ix_i}{\sum w_i} \f$
88    ///
89    double mean(void) const;
90
91    ///
92    /// @brief Weighted version of number of data points.
93    ///
94    /// If all
95    /// weights are equal, the unweighted version is identical to the
96    /// non-weighted version. Adding a data point with zero weight
97    /// does not change n(). The calculated value is always smaller
98    /// than the actual number of data points added to object.
99    ///
100    /// @return \f$ \frac{\left(\sum w_i\right)^2}{\sum w_i^2} \f$
101    ///
102    double n(void) const;
103
104    ///
105    /// @brief Rescale object.
106    ///
107    /// Each data point is rescaled as \f$ x = a * x \f$
108    ///
109    void rescale(double a);
110
111    ///
112    /// @brief Reset everything to zero.
113    ///
114    void reset(void);
115
116    ///
117    /// @brief The standard deviation is defined as the square root of
118    /// the variance().
119    ///
120    /// @return The standard deviation, root of the variance().
121    ///
122    double std(void) const;
123
124    ///
125    /// @brief Calculates standard deviation of the mean().
126    ///
127    /// Variance from the
128    /// weights are here neglected. This is true when the weight is
129    /// known before the measurement. In case this is not a good
130    /// approximation, use bootstrapping to estimate the error.
131    ///
132    /// @return \f$ \frac{\sum w^2}{\left(\sum w\right)^3}\sum w(x-m)^2 \f$
133    /// where \f$ m \f$ is the mean()
134    ///
135    double standard_error(void) const;
136
137    ///
138    /// Calculating the sum of weights 
139    ///
140    /// @return \f$ \sum w_i \f$
141    ///
142    double sum_w(void) const;
143
144    ///
145    /// @return \f$ \sum w_i^2 \f$
146    ///
147    double sum_ww(void) const;
148
149    ///
150    /// \f$ \sum w_ix_i \f$
151    ///
152    /// @return weighted sum of x
153    ///
154    double sum_wx(void) const;
155
156    ///
157    ///  @return \f$ \sum w_i x_i^2 \f$
158    ///
159    double sum_wxx(void) const;
160
161    ///
162    /// @return \f$ \sum_i w_i (x_i-m)^2\f$
163    ///
164    double sum_xx_centered(void) const;
165
166    /**
167       The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
168       }{\sum w_i} \f$, where \a m is the known mean.
169       
170       @return Variance when the mean is known to be \a m.
171    */
172    double variance(const double m) const;
173
174    /**
175       The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
176       }{\sum w_i} \f$, where \a m is the mean(). Here the weight are
177       interpreted as probability weights. For analytical weights the
178       variance has no meaning as each data point has its own
179       variance.
180       
181       @return The variance.
182    */
183    double variance(void) const;
184
185    /**
186       operator to add an AveragerWeighted
187
188       \since New in yat 0.6
189    */
190    const AveragerWeighted& operator+=(const AveragerWeighted&);
191   
192  private:
193    ///
194    ///  @return \f$ \sum w_i^2x_i \f$
195    ///
196    double sum_wwx(void) const;
197
198    ///
199    ///  @return \f$ \sum w_i^2x_i^2 \f$
200    ///
201    double sum_wwxx(void) const;
202   
203    double mean_; // wx/w
204    double m2_; // wxx - m*m*w
205    double w_;
206    double ww_;
207    double wwxx_;
208    double wxx_;
209
210    // double wwx; // m2_ + m*m*w
211    // double wx; // w*mean_
212  };
213
214  /**
215     \brief adding a range of values to AveragerWeighted \a a
216
217     If InputIterator is non-weighted unitary weights are used.
218
219     \relates AveragerWeighted
220   */
221  template <typename InputIterator>
222  void add(AveragerWeighted& a, InputIterator first, InputIterator last)
223  {
224    BOOST_CONCEPT_ASSERT((boost::InputIterator<InputIterator>));
225    utility::iterator_traits<InputIterator> traits;
226    for ( ; first != last; ++first)
227      a.add(traits.data(first), traits.weight(first));
228  }
229
230  /**
231     \brief add values from two ranges to AveragerWeighted \a a
232
233     Add data from range [first1, last1) with their corresponding
234     weight in range [first2, first2 + distance(first, last) ).
235
236     Requirement: InputIterator1 and InputIterator2 are unweighted
237     input iterators.
238
239     \relates AveragerWeighted
240   */
241  template <typename InputIterator1, typename InputIterator2>
242  void add(AveragerWeighted& a, InputIterator1 first1, InputIterator1 last1, 
243           InputIterator2 first2)
244  {
245    BOOST_CONCEPT_ASSERT((boost::InputIterator<InputIterator1>));
246    BOOST_CONCEPT_ASSERT((boost::InputIterator<InputIterator2>));
247    utility::check_iterator_is_unweighted(first1);
248    utility::check_iterator_is_unweighted(first2);
249    for ( ; first1 != last1; ++first1, ++first2)
250      a.add(*first1, *first2);
251  }
252
253}}} // of namespace statistics, yat, and theplu
254
255#endif
Note: See TracBrowser for help on using the repository browser.