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

Last change on this file since 914 was 914, checked in by Peter, 14 years ago

refactoring add function in Averagers, refs #246

  • 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 914 2007-09-29 02:18:52Z peter $
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/trac/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/IteratorTraits.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    /// Adding each value in an array \a x and corresponding value in
86    /// weight array \a w.
87    ///
88    /// The requirements for the types T1 and T2 of the arrays \a x
89    /// and \a w are: operator[] returning an element and function
90    /// size() returning the number of elements.
91    ///
92    template <typename T1, typename T2>
93    void add_values(const T1& x, const T2& w);
94
95    ///
96    /// @brief Calculate the weighted mean
97    ///
98    /// @return \f$ \frac{\sum w_ix_i}{\sum w_i} \f$
99    ///
100    double mean(void) const;
101
102    ///
103    /// @brief Weighted version of number of data points.
104    ///
105    /// If all
106    /// weights are equal, the unweighted version is identical to the
107    /// non-weighted version. Adding a data point with zero weight
108    /// does not change n(). The calculated value is always smaller
109    /// than the actual number of data points added to object.
110    ///
111    /// @return \f$ \frac{\left(\sum w_i\right)^2}{\sum w_i^2} \f$
112    ///
113    double n(void) const;
114
115    ///
116    /// @brief Rescale object.
117    ///
118    /// Each data point is rescaled as \f$ x = a * x \f$
119    ///
120    void rescale(double a);
121
122    ///
123    /// @brief Reset everything to zero.
124    ///
125    void reset(void);
126
127    ///
128    /// @brief The standard deviation is defined as the square root of
129    /// the variance().
130    ///
131    /// @return The standard deviation, root of the variance().
132    ///
133    double std(void) const;
134
135    ///
136    /// @brief Calculates standard deviation of the mean().
137    ///
138    /// Variance from the
139    /// weights are here neglected. This is true when the weight is
140    /// known before the measurement. In case this is not a good
141    /// approximation, use bootstrapping to estimate the error.
142    ///
143    /// @return \f$ \frac{\sum w^2}{\left(\sum w\right)^3}\sum w(x-m)^2 \f$
144    /// where \f$ m \f$ is the mean()
145    ///
146    double standard_error(void) const;
147
148    ///
149    /// Calculating the sum of weights 
150    ///
151    /// @return \f$ \sum w_i \f$
152    ///
153    double sum_w(void) const;
154
155    ///
156    /// @return \f$ \sum w_i^2 \f$
157    ///
158    double sum_ww(void) const;
159
160    ///
161    /// \f$ \sum w_ix_i \f$
162    ///
163    /// @return weighted sum of x
164    ///
165    double sum_wx(void) const;
166
167    ///
168    ///  @return \f$ \sum w_i x_i^2 \f$
169    ///
170    double sum_wxx(void) const;
171
172    ///
173    /// @return \f$ \sum_i w_i (x_i-m)^2\f$
174    ///
175    double sum_xx_centered(void) const;
176
177    /**
178       The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
179       }{\sum w_i} \f$, where \a m is the known mean.
180       
181       @return Variance when the mean is known to be \a m.
182    */
183    double variance(const double m) const;
184
185    /**
186       The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
187       }{\sum w_i} \f$, where \a m is the mean(). Here the weight are
188       interpreted as probability weights. For analytical weights the
189       variance has no meaning as each data point has its own
190       variance.
191       
192       @return The variance.
193    */
194    double variance(void) const;
195
196
197  private:
198    ///
199    ///  @return \f$ \sum w_i^2x_i \f$
200    ///
201    double sum_wwx(void) const;
202
203    ///
204    ///  @return \f$ \sum w_i^2x_i^2 \f$
205    ///
206    double sum_wwxx(void) const;
207   
208    const Averager& wx(void) const;
209    const Averager& w(void) const;
210
211    ///
212    /// operator to add a AveragerWeighted
213    ///
214    const AveragerWeighted& operator+=(const AveragerWeighted&);
215   
216    Averager w_;
217    Averager wx_;
218    double wwx_;
219    double wxx_;
220  };
221
222  /**
223     \brief adding a ranges of values to Averager \a a
224   */
225  template <typename Iter>
226  void add(AveragerWeighted& a, Iter first, Iter last)
227  {
228    for ( ; first != last; ++first)
229      a.add(utility::iterator_traits_data(first));
230  }
231
232  // Template implementations
233  template <typename T1, typename T2>
234  void AveragerWeighted::add_values(const T1& x, const T2& w)
235  {
236    for (size_t i=0; i<x.size(); i++) 
237      add(x[i],w[i]);
238  }
239
240///
241/// The AveragerWeighted output operator
242///
243///std::ostream& operator<<(std::ostream& s,const AveragerWeighted&);
244
245}}} // of namespace statistics, yat, and theplu
246
247#endif
Note: See TracBrowser for help on using the repository browser.