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

Last change on this file since 1437 was 1437, checked in by Peter, 13 years ago

merge patch release 0.4.2 to trunk. Delta 0.4.2-0.4.1

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.7 KB
Line 
1#ifndef _theplu_yat_statistics_averagerweighted_
2#define _theplu_yat_statistics_averagerweighted_
3
4// $Id: AveragerWeighted.h 1437 2008-08-25 17:55:00Z 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 Jari Häkkinen, Peter Johansson
10  Copyright (C) 2008 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 2 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 this program; if not, write to the Free Software
26  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
27  02111-1307, USA.
28*/
29
30#include "Averager.h"
31#include "yat/utility/iterator_traits.h"
32
33#include <cmath>
34
35namespace theplu{
36namespace yat{
37namespace statistics{
38
39  ///
40  /// @brief Class to calulate averages with weights.
41  ///
42  /// There are several different reasons why a statistical analysis
43  /// needs to adjust for weighting. In the litterature reasons are
44  /// mainly divided into two kinds of weights - probablity weights
45  /// and analytical weights. 1) Analytical weights are appropriate
46  /// for scientific experiments where some measurements are known to
47  /// be more precise than others. The larger weight a measurement has
48  /// the more precise is is assumed to be, or more formally the
49  /// weight is proportional to the reciprocal variance
50  /// \f$ \sigma_i^2 = \frac{\sigma^2}{w_i} \f$. 2) Probability weights
51  /// are used for the situation when calculating averages over a
52  /// distribution \f$ f \f$ , but sampling from a distribution \f$ f' \f$.
53  /// Compensating for this discrepancy averages of observables
54  /// are taken to be \f$ \sum \frac{f}{f'}X \f$ For further discussion:
55  /// see \ref weighted_statistics
56  ///
57  /// If nothing else stated, each function fulfills the
58  /// following:<br> <ul><li>Setting a weight to zero corresponds to
59  /// removing the data point from the dataset.</li><li> Setting all
60  /// weights to unity, the yields the same result as from
61  /// corresponding function in Averager.</li><li> Rescaling weights
62  /// does not change the performance of the object.</li></ul>
63  ///
64  /// @see Averager AveragerPair AveragerPairWeighted
65  ///
66  class AveragerWeighted
67  {
68  public:
69
70    ///
71    /// @brief The default constructor
72    ///
73    AveragerWeighted(void);
74
75    ///
76    /// @brief The copy constructor
77    ///
78    AveragerWeighted(const AveragerWeighted&);
79
80    ///
81    /// Adding a data point \a d, with weight \a w (default is 1)
82    ///
83    void  add(const double d,const double w=1);
84
85    ///
86    /// @brief Calculate the weighted mean
87    ///
88    /// @return \f$ \frac{\sum w_ix_i}{\sum w_i} \f$
89    ///
90    double mean(void) const;
91
92    ///
93    /// @brief Weighted version of number of data points.
94    ///
95    /// If all
96    /// weights are equal, the unweighted version is identical to the
97    /// non-weighted version. Adding a data point with zero weight
98    /// does not change n(). The calculated value is always smaller
99    /// than the actual number of data points added to object.
100    ///
101    /// @return \f$ \frac{\left(\sum w_i\right)^2}{\sum w_i^2} \f$
102    ///
103    double n(void) const;
104
105    ///
106    /// @brief Rescale object.
107    ///
108    /// Each data point is rescaled as \f$ x = a * x \f$
109    ///
110    void rescale(double a);
111
112    ///
113    /// @brief Reset everything to zero.
114    ///
115    void reset(void);
116
117    ///
118    /// @brief The standard deviation is defined as the square root of
119    /// the variance().
120    ///
121    /// @return The standard deviation, root of the variance().
122    ///
123    double std(void) const;
124
125    ///
126    /// @brief Calculates standard deviation of the mean().
127    ///
128    /// Variance from the
129    /// weights are here neglected. This is true when the weight is
130    /// known before the measurement. In case this is not a good
131    /// approximation, use bootstrapping to estimate the error.
132    ///
133    /// @return \f$ \frac{\sum w^2}{\left(\sum w\right)^3}\sum w(x-m)^2 \f$
134    /// where \f$ m \f$ is the mean()
135    ///
136    double standard_error(void) const;
137
138    ///
139    /// Calculating the sum of weights 
140    ///
141    /// @return \f$ \sum w_i \f$
142    ///
143    double sum_w(void) const;
144
145    ///
146    /// @return \f$ \sum w_i^2 \f$
147    ///
148    double sum_ww(void) const;
149
150    ///
151    /// \f$ \sum w_ix_i \f$
152    ///
153    /// @return weighted sum of x
154    ///
155    double sum_wx(void) const;
156
157    ///
158    ///  @return \f$ \sum w_i x_i^2 \f$
159    ///
160    double sum_wxx(void) const;
161
162    ///
163    /// @return \f$ \sum_i w_i (x_i-m)^2\f$
164    ///
165    double sum_xx_centered(void) const;
166
167    /**
168       The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
169       }{\sum w_i} \f$, where \a m is the known mean.
170       
171       @return Variance when the mean is known to be \a m.
172    */
173    double variance(const double m) const;
174
175    /**
176       The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
177       }{\sum w_i} \f$, where \a m is the mean(). Here the weight are
178       interpreted as probability weights. For analytical weights the
179       variance has no meaning as each data point has its own
180       variance.
181       
182       @return The variance.
183    */
184    double variance(void) const;
185
186
187  private:
188    ///
189    ///  @return \f$ \sum w_i^2x_i \f$
190    ///
191    double sum_wwx(void) const;
192
193    ///
194    ///  @return \f$ \sum w_i^2x_i^2 \f$
195    ///
196    double sum_wwxx(void) const;
197   
198    const Averager& wx(void) const;
199    const Averager& w(void) const;
200
201    ///
202    /// operator to add a AveragerWeighted
203    ///
204    const AveragerWeighted& operator+=(const AveragerWeighted&);
205   
206    Averager w_;
207    Averager wx_;
208    double wwx_;
209    double wxx_;
210  };
211
212  /**
213     \brief adding a range of values to AveragerWeighted \a a
214
215     If iterator is non-weighted unitary weights are used.
216   */
217  template <typename Iter>
218  void add(AveragerWeighted& a, Iter first, Iter last)
219  {
220    for ( ; first != last; ++first)
221      a.add(utility::iterator_traits<Iter>().data(first),
222            utility::iterator_traits<Iter>().weight(first));
223  }
224
225  /**
226     \brief add values from two ranges to AveragerWeighted \a a
227
228     Add data from range [first1, last1) with their corresponding
229     weight in range [first2, first2 + distance(first, last) ).
230
231     Requirement: Iter1 and Iter2 are unweighted iterators.
232   */
233  template <typename Iter1, typename Iter2>
234  void add(AveragerWeighted& a, Iter1 first1, Iter1 last1, Iter2 first2)
235  {
236    utility::check_iterator_is_unweighted(first1);
237    utility::check_iterator_is_unweighted(first2);
238    for ( ; first1 != last1; ++first1, ++first2)
239      a.add(*first1, *first2);
240  }
241
242}}} // of namespace statistics, yat, and theplu
243
244#endif
Note: See TracBrowser for help on using the repository browser.