# source:trunk/lib/statistics/AveragerWeighted.h@494

Last change on this file since 494 was 494, checked in by Peter, 17 years ago

• Property svn:eol-style set to native
• Property svn:keywords set to Author Date Id Revision
File size: 5.7 KB
Line
1// $Id: AveragerWeighted.h 494 2006-01-10 13:44:14Z peter$
2
3#ifndef _theplu_statistics_averager_weighted_
4#define _theplu_statistics_averager_weighted_
5
6#include <c++_tools/statistics/Averager.h>
7
8#include <cmath>
9//#include <ostream>
10
11namespace theplu{
12  class gslapi::vector;
13
14namespace statistics{
15
16  ///
17  /// @brief Class to calulate averages with weights.
18  ///
19  /// There are several different reasons why a statistical analysis
20  /// needs to adjust for weighting. In the litterature reasons are
21  /// mainly divided into two kinds of weights - probablity weights
22  /// and analytical weights. 1) Analytical weights are appropriate
23  /// for scientific experiments where some measurements are known to
24  /// be more precise than others. The larger weight a measurement has
25  /// the more precise is is assumed to be, or more formally the
26  /// weight is proportional to the reciprocal variance
27  /// \f$\sigma_i^2 = \frac{\sigma^2}{w_i} \f$. 2) Probablity weights
28  /// are used for the situation when calculating averages over a
29  /// distribution \f$f \f$ , but sampling from a distribution \f$f' 30 /// \f$. Compensating for this discrepancy averages of observables
31  /// are taken to be \f$\sum \frac{f}{f'}X \f$ For further discussion:
32  /// <a href="Statistics/index.html">Weighted Statistics document</a><br>
33  ///
34  /// If nothing else stated, each function fulfills the
35  /// following:<br> <ul><li>Setting a weight to zero corresponds to
36  /// removing the data point from the dataset.</li><li> Setting all
37  /// weights to unity, the yields the same result as from
38  /// corresponding function in Averager.</li><li> Rescaling weights
39  /// does not change the performance of the object.</li></ul>
40  ///
41  /// @see Averager AveragerPair AveragerPairWeighted
42  ///
43  class AveragerWeighted
44  {
45  public:
46
47    ///
48    /// Default constructor
49    ///
50    inline AveragerWeighted(void)
51      : w_(Averager()), wx_(Averager()), wwx_(0), wxx_(0) {}
52
53    ///
54    /// Copy constructor
55    ///
56    inline AveragerWeighted(const AveragerWeighted& a)
57      : w_(Averager(a.sum_w(),a.sum_ww(),1)),
58        wx_(Averager(a.sum_wx(),a.sum_wwxx(),1)), wwx_(a.sum_wwx()),
59        wxx_(a.sum_wxx()) {}
60
61    ///
62    /// adding a data point d, with weight w (default is 1)
63    ///
64    inline void add(const double d,const double w=1)
66
67    ///
68    /// Adding each value in vector \a x and corresponding value in
69    /// weight vector \a w.
70    ///
71    void add(const gslapi::vector& x, const gslapi::vector& w);
72
73    ///
74    /// Calculating the weighted mean
75    ///
76    /// @return \f$\frac{\sum w_ix_i}{\sum w_i} \f$
77    ///
78    inline double mean(void) const { return sum_w() ?
79                                       sum_wx()/sum_w() : 0; }
80
81    ///
82    /// @return \f$\frac{\left(\sum w_i\right)^2}{\sum w_i^2} \f$
83    ///
84    inline double n(void) const { return sum_w()*sum_w()/sum_ww(); }
85
86    ///
87    /// rescale object, i.e. each data point is rescaled
88    /// \f$x = a * x \f$
89    ///
90    inline void rescale(double a) { wx_.rescale(a); wwx_*=a; wxx_*=a*a; }
91
92    ///
93    /// resets everything to zero
94    ///
95    inline void reset(void) { wx_.reset(); w_.reset(); wwx_=0; wxx_=0; }
96
97    ///
98    /// The standard deviation is defined as the square root of the
99    /// variance().
100    ///
101    /// @return The standard deviation, root of the variance().
102    ///
103    inline double std(void) const { return sqrt(variance()); }
104
105    ///
106    /// Calculates standard deviation of the mean(). Variance from the
107    /// weights are here neglected. This is true when the weight is
108    /// known before the measurement. In case this is not a good
109    /// approximation, use bootstrapping to estimate the error.
110    ///
111    /// @return \f$\frac{\sum w^2}{\left(\sum w\right)^3}\sum w(x-m)^2 \f$
112    /// where \f$m \f$ is the mean()
113    ///
114    inline double standard_error(void)  const
115    { return sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) *
116                  sum_xx_centered()); }
117
118    ///
119    /// Calculating the sum of weights
120    ///
121    /// @return \f$\sum w_i \f$
122    ///
123    inline double sum_w(void) const
124    { return w_.sum_x(); }
125
126    ///
127    /// @return \f$\sum w_i^2 \f$
128    ///
129    inline double sum_ww(void)  const
130    { return w_.sum_xx(); }
131
132    ///
133    /// \f$\sum w_ix_i \f$
134    ///
135    /// @return weighted sum of x
136    ///
137    inline double sum_wx(void)  const
138    { return wx_.sum_x(); }
139
140    ///
141    /// @return \f$\sum_i w_i (x_i-m)^2\f$
142    ///
143    inline double sum_xx_centered(void) const
144    { return sum_wxx() - mean()*mean()*sum_w(); }
145
146    ///
147    /// The variance is calculated as \f$\frac{\sum w_i (x_i - m)^2 148 /// }{\sum w_i} \f$, where \a m is the known mean.
149    ///
150    /// @return Variance when the mean is known to be \a m.
151    ///
152    inline double variance(const double m) const
153    { return (sum_wxx()-2*m*sum_wx())/sum_w()+m*m; }
154
155    ///
156    /// The variance is calculated as \f$\frac{\sum w_i (x_i - m)^2 157 /// }{\sum w_i} \f$, where \a m is the mean(). Here the weight are
158    /// interpreted as probability weights. For analytical weights the
159    /// variance has no meaning as each data point has its own
160    /// variance.
161    ///
162    /// @return The variance.
163    ///
164    inline double variance(void) const
165    { return sum_xx_centered()/sum_w(); }
166
167
168  private:
169    ///
170    ///  @return \f$\sum w_i^2x_i^2 \f$
171    ///
172    inline double sum_wwxx(void)  const
173    { return wx_.sum_xx(); }
174
175    ///
176    ///  @return \f$\sum w_i^2x_i \f$
177    ///
178    inline double sum_wwx(void) const
179    { return wwx_; }
180
181    ///
182    ///  @return \f$\sum w_i x_i^2 \f$
183    ///
184    inline double sum_wxx(void) const { return wxx_; }
185
186    ///
187    /// operator to add a AveragerWeighted
188    ///
189    AveragerWeighted operator+=(const AveragerWeighted&);
190
191    Averager w_;
192    Averager wx_;
193    double wwx_;
194    double wxx_;
195
196    inline Averager wx(void) const {return wx_;}
197    inline Averager w(void) const {return w_;}
198
199
200  };
201
202///
203/// The AveragerWeighted output operator
204///
205///std::ostream& operator<<(std::ostream& s,const AveragerWeighted&);
206
207}} // of namespace statistics and namespace theplu
208
209#endif
Note: See TracBrowser for help on using the repository browser.