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

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

changing URL to http://trac.thep.lu.se/trac/yat

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