source: trunk/yat/statistics/AveragerPairWeighted.h @ 1122

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

refs #335 - Changed so Averager classes are consistently returning NaN when Averager is empty or for some other reason the estimation ends up with things like zero by zero division. Previously zero was returned from some functions and Nan from some functions. I did not change anythuing in NCC.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.8 KB
Line 
1#ifndef _theplu_yat_statistics_averagerpairweighted_
2#define _theplu_yat_statistics_averagerpairweighted_
3
4// $Id: AveragerPairWeighted.h 1122 2008-02-22 17:01:15Z peter $
5
6/*
7  Copyright (C) 2005 Markus Ringnér, Peter Johansson
8  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
9  Copyright (C) 2007 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 "AveragerWeighted.h"
30
31#include "yat/utility/iterator_traits.h"
32#include "yat/utility/yat_assert.h"
33
34#include <cmath>
35#include <stdexcept>
36
37namespace theplu{
38namespace yat{
39namespace statistics{
40  ///
41  /// @brief Class for taking care of mean and covariance of two variables in
42  /// a weighted manner.
43  ///
44  /// <a href="Statistics/index.html">Weighted Statistics document</a>
45  ///
46  /// If nothing else stated, each function fulfills the
47  /// following:<br> <ul><li>Setting a weight to zero corresponds to
48  /// removing the data point from the dataset.</li><li> Setting all
49  /// weights to unity, the yields the same result as from
50  /// corresponding function in AveragerPair.</li><li> Rescaling weights
51  /// does not change the performance of the object.</li></ul>
52  ///
53  /// @see Averager AveragerWeighted AveragerPair
54  ///
55  class AveragerPairWeighted
56  {
57  public:
58
59    ///
60    /// @brief The default constructor
61    ///
62    AveragerPairWeighted(void);
63
64    ///
65    /// Adding a pair of data points with value \a x and \a y, and
66    /// their weights. If either of the weights are zero the addition
67    /// is ignored
68    ///
69    void  add(const double x, const double y, 
70              const double wx, const double wy);
71
72    ///
73    /// @brief Pearson correlation coefficient.
74    ///
75    /// @return \f$ \frac{\sum w_xw_y (x-m_x)(y-m_y)}{\sqrt{\sum
76    /// w_xw_y (y-m_y)^2\sum w_xw_y (y-m_y)^2}} \f$ where m is
77    /// calculated as \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
78    ///
79    double correlation(void) const;
80 
81    ///
82    /// \f$ \frac{\sum w_xw_y (x-m_x)(y-m_y)}{\sum w_xw_y} \f$ where m
83    /// is calculated as \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
84    ///
85    double covariance(void) const;
86
87    /**
88       @return \f$ \frac{\sum w_xw_y(x-y)^2}{\sum w_xw_y} \f$
89    */
90    double msd(void) const;
91
92    /**
93       @return \f$ \frac{\left(\sum w_x w_y\right)^2}{\sum w_x^2w_y^2} \f$
94    */ 
95    double n(void) const;
96
97    ///
98    /// @brief Reset everything to zero
99    ///
100    void reset(void);
101
102    ///
103    /// @return \f$ \sum w_xw_y \f$
104    ///
105    double sum_w(void) const;
106
107    ///
108    /// @return \f$ \sum w_xw_yxy \f$
109    ///
110    double sum_xy(void) const;
111
112    ///
113    /// @return \f$ \sum w_xw_y (x-m_x)(y-m_y) \f$ where m is calculated as
114    /// \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
115    ///
116    double sum_xy_centered(void) const;
117
118    ///
119    /// @note the weights are calculated as \f$ w = w_x * w_y \f$
120    ///
121    /// @return AveragerWeighted for x
122    ///
123    const AveragerWeighted& x_averager(void) const;
124
125    ///
126    /// @note the weights are calculated as \f$ w = w_x * w_y \f$
127    ///
128    /// @return AveragerWeighted for y
129    ///
130    const AveragerWeighted& y_averager(void) const;
131
132  private:
133    AveragerWeighted x_; // weighted averager with w = w_x*w_y
134    AveragerWeighted y_; // weighted averager with w = w_x*w_y
135    double wxy_;
136    double w_;
137
138  };
139
140  /**
141     \brief adding a ranges of values to AveragerPairWeighted \a ap
142   */
143  template <class Iter1, class Iter2>
144  void add(AveragerPairWeighted& ap, Iter1 first1, Iter1 last1, Iter2 first2)
145  {
146    for ( ; first1 != last1; ++first1, ++first2) {
147      ap.add(utility::iterator_traits<Iter1>().data(first1), 
148             utility::iterator_traits<Iter2>().data(first2),
149             utility::iterator_traits<Iter1>().weight(first1), 
150             utility::iterator_traits<Iter2>().weight(first2));
151    }
152  }
153
154
155  /**
156     \brief adding four ranges of values to AveragerPairWeighted \a ap
157   */
158  template <class Iter1, class Iter2, class Iter3, class Iter4>
159  void add(AveragerPairWeighted& ap, Iter1 x, Iter1 xlast, Iter2 y, Iter3 wx, 
160           Iter4 wy)
161  {
162    utility::check_iterator_is_unweighted(x);
163    utility::check_iterator_is_unweighted(y);
164    utility::check_iterator_is_unweighted(wx);
165    utility::check_iterator_is_unweighted(wy);
166    while (x!=xlast) {
167      ap.add(*x, *y, *wx, *wy);
168      ++x;
169      ++y;
170      ++wx;
171      ++wy;
172    }
173  }
174
175
176}}} // of namespace statistics, yat, and theplu
177
178#endif
Note: See TracBrowser for help on using the repository browser.