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

Last change on this file since 1000 was 1000, checked in by Jari Häkkinen, 14 years ago

trac moved to new location.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.6 KB
Line 
1#ifndef _theplu_yat_statistics_averagerweighted_
2#define _theplu_yat_statistics_averagerweighted_
3
4// $Id: AveragerWeighted.h 1000 2007-12-23 20:09:15Z jari $
5
6/*
7  Copyright (C) 2004 Cecilia Ritz, Jari Häkkinen, Peter Johansson
8  Copyright (C) 2005, 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
9  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
10
11  This file is part of the yat library, http://trac.thep.lu.se/yat
12
13  The yat library is free software; you can redistribute it and/or
14  modify it under the terms of the GNU General Public License as
15  published by the Free Software Foundation; either version 2 of the
16  License, or (at your option) any later version.
17
18  The yat library is distributed in the hope that it will be useful,
19  but WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21  General Public License for more details.
22
23  You should have received a copy of the GNU General Public License
24  along with this program; if not, write to the Free Software
25  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
26  02111-1307, USA.
27*/
28
29#include "Averager.h"
30#include "yat/utility/iterator_traits.h"
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  /// <a href="Statistics/index.html">Weighted Statistics document</a><br>
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  private:
187    ///
188    ///  @return \f$ \sum w_i^2x_i \f$
189    ///
190    double sum_wwx(void) const;
191
192    ///
193    ///  @return \f$ \sum w_i^2x_i^2 \f$
194    ///
195    double sum_wwxx(void) const;
196   
197    const Averager& wx(void) const;
198    const Averager& w(void) const;
199
200    ///
201    /// operator to add a AveragerWeighted
202    ///
203    const AveragerWeighted& operator+=(const AveragerWeighted&);
204   
205    Averager w_;
206    Averager wx_;
207    double wwx_;
208    double wxx_;
209  };
210
211  /**
212     \brief adding a range of values to AveragerWeighted \a a
213   */
214  template <typename Iter>
215  void add(AveragerWeighted& a, Iter first, Iter last)
216  {
217    for ( ; first != last; ++first)
218      a.add(utility::iterator_traits_data(first),
219            utility::iterator_traits_weight(first));
220  }
221
222  /**
223     \brief add values from two ranges to AveragerWeighted \a a
224
225     Add data from range [first1, last1) with their corresponding
226     weight in range [first2, first2 + distance(first, last) ).
227
228     Requirement: Iter1 and Iter2 are unweighted iterators.
229   */
230  template <typename Iter1, typename Iter2>
231  void add(AveragerWeighted& a, Iter1 first1, Iter1 last1, Iter2 first2)
232  {
233    utility::check_iterator_is_unweighted(first1);
234    utility::check_iterator_is_unweighted(first2);
235    for ( ; first1 != last1; ++first1, ++first2)
236      a.add(*first1, *first2);
237  }
238
239}}} // of namespace statistics, yat, and theplu
240
241#endif
Note: See TracBrowser for help on using the repository browser.