source: trunk/lib/statistics/AveragerWeighted.h @ 492

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

updated Weighted Statistics document and reimplemented few stuff. NOTE, this revision may not compile. Cannot test compilation on Le Mac.

  • 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 492 2006-01-09 09:38:12Z 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)
65    { if(w==0) return; w_.add(w); wx_.add(w*d); wwx_+=w*w*d; wxx_+=w*d*d; }
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.