yat/statistics/AveragerPairWeighted.h

Code
Comments
Other
Rev Date Author Line
683 11 Oct 06 jari 1 #ifndef _theplu_yat_statistics_averagerpairweighted_
683 11 Oct 06 jari 2 #define _theplu_yat_statistics_averagerpairweighted_
474 22 Dec 05 peter 3
616 31 Aug 06 jari 4 // $Id$
616 31 Aug 06 jari 5
675 10 Oct 06 jari 6 /*
2119 12 Dec 09 peter 7   Copyright (C) 2005 Peter Johansson, Markus Ringnér
4359 23 Aug 23 peter 8   Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
4359 23 Aug 23 peter 9   Copyright (C) 2007 Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 10   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 11   Copyright (C) 2009, 2010, 2011, 2016 Peter Johansson
474 22 Dec 05 peter 12
1437 25 Aug 08 peter 13   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 14
675 10 Oct 06 jari 15   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 16   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 17   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 18   License, or (at your option) any later version.
675 10 Oct 06 jari 19
675 10 Oct 06 jari 20   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 21   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 22   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 23   General Public License for more details.
675 10 Oct 06 jari 24
675 10 Oct 06 jari 25   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 26   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 27 */
675 10 Oct 06 jari 28
680 11 Oct 06 jari 29 #include "AveragerWeighted.h"
675 10 Oct 06 jari 30
937 05 Oct 07 peter 31 #include "yat/utility/iterator_traits.h"
916 30 Sep 07 peter 32 #include "yat/utility/yat_assert.h"
914 29 Sep 07 peter 33
2202 21 Feb 10 peter 34 #include <boost/concept_check.hpp>
2202 21 Feb 10 peter 35
474 22 Dec 05 peter 36 #include <cmath>
936 05 Oct 07 peter 37 #include <stdexcept>
474 22 Dec 05 peter 38
474 22 Dec 05 peter 39 namespace theplu{
680 11 Oct 06 jari 40 namespace yat{
474 22 Dec 05 peter 41 namespace statistics{
474 22 Dec 05 peter 42   ///
767 22 Feb 07 peter 43   /// @brief Class for taking care of mean and covariance of two variables in
474 22 Dec 05 peter 44   /// a weighted manner.
474 22 Dec 05 peter 45   ///
1185 28 Feb 08 peter 46   /// \see \ref weighted_statistics
489 04 Jan 06 peter 47   ///
489 04 Jan 06 peter 48   /// If nothing else stated, each function fulfills the
489 04 Jan 06 peter 49   /// following:<br> <ul><li>Setting a weight to zero corresponds to
489 04 Jan 06 peter 50   /// removing the data point from the dataset.</li><li> Setting all
489 04 Jan 06 peter 51   /// weights to unity, the yields the same result as from
489 04 Jan 06 peter 52   /// corresponding function in AveragerPair.</li><li> Rescaling weights
489 04 Jan 06 peter 53   /// does not change the performance of the object.</li></ul>
489 04 Jan 06 peter 54   ///
490 04 Jan 06 peter 55   /// @see Averager AveragerWeighted AveragerPair
489 04 Jan 06 peter 56   ///
474 22 Dec 05 peter 57   class AveragerPairWeighted
474 22 Dec 05 peter 58   {
474 22 Dec 05 peter 59   public:
474 22 Dec 05 peter 60
703 18 Dec 06 jari 61     ///
703 18 Dec 06 jari 62     /// @brief The default constructor
703 18 Dec 06 jari 63     ///
703 18 Dec 06 jari 64     AveragerPairWeighted(void);
476 22 Dec 05 markus 65
474 22 Dec 05 peter 66     ///
474 22 Dec 05 peter 67     /// Adding a pair of data points with value \a x and \a y, and
489 04 Jan 06 peter 68     /// their weights. If either of the weights are zero the addition
489 04 Jan 06 peter 69     /// is ignored
474 22 Dec 05 peter 70     ///
4200 19 Aug 22 peter 71     void  add(const double x, const double y,
474 22 Dec 05 peter 72               const double wx, const double wy);
474 22 Dec 05 peter 73
474 22 Dec 05 peter 74     ///
1186 28 Feb 08 peter 75     /// @brief %Pearson correlation coefficient.
474 22 Dec 05 peter 76     ///
489 04 Jan 06 peter 77     /// @return \f$ \frac{\sum w_xw_y (x-m_x)(y-m_y)}{\sqrt{\sum
489 04 Jan 06 peter 78     /// w_xw_y (y-m_y)^2\sum w_xw_y (y-m_y)^2}} \f$ where m is
489 04 Jan 06 peter 79     /// calculated as \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
489 04 Jan 06 peter 80     ///
718 26 Dec 06 jari 81     double correlation(void) const;
4200 19 Aug 22 peter 82
474 22 Dec 05 peter 83     ///
489 04 Jan 06 peter 84     /// \f$ \frac{\sum w_xw_y (x-m_x)(y-m_y)}{\sum w_xw_y} \f$ where m
489 04 Jan 06 peter 85     /// is calculated as \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
474 22 Dec 05 peter 86     ///
718 26 Dec 06 jari 87     double covariance(void) const;
489 04 Jan 06 peter 88
772 25 Feb 07 peter 89     /**
772 25 Feb 07 peter 90        @return \f$ \frac{\sum w_xw_y(x-y)^2}{\sum w_xw_y} \f$
772 25 Feb 07 peter 91     */
772 25 Feb 07 peter 92     double msd(void) const;
772 25 Feb 07 peter 93
899 26 Sep 07 peter 94     /**
899 26 Sep 07 peter 95        @return \f$ \frac{\left(\sum w_x w_y\right)^2}{\sum w_x^2w_y^2} \f$
4200 19 Aug 22 peter 96     */
899 26 Sep 07 peter 97     double n(void) const;
899 26 Sep 07 peter 98
489 04 Jan 06 peter 99     ///
703 18 Dec 06 jari 100     /// @brief Reset everything to zero
489 04 Jan 06 peter 101     ///
703 18 Dec 06 jari 102     void reset(void);
474 22 Dec 05 peter 103
890 25 Sep 07 markus 104     ///
489 04 Jan 06 peter 105     /// @return \f$ \sum w_xw_y \f$
489 04 Jan 06 peter 106     ///
718 26 Dec 06 jari 107     double sum_w(void) const;
489 04 Jan 06 peter 108
489 04 Jan 06 peter 109     ///
489 04 Jan 06 peter 110     /// @return \f$ \sum w_xw_yxy \f$
489 04 Jan 06 peter 111     ///
718 26 Dec 06 jari 112     double sum_xy(void) const;
489 04 Jan 06 peter 113
489 04 Jan 06 peter 114     ///
489 04 Jan 06 peter 115     /// @return \f$ \sum w_xw_y (x-m_x)(y-m_y) \f$ where m is calculated as
489 04 Jan 06 peter 116     /// \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
489 04 Jan 06 peter 117     ///
718 26 Dec 06 jari 118     double sum_xy_centered(void) const;
489 04 Jan 06 peter 119
489 04 Jan 06 peter 120     ///
489 04 Jan 06 peter 121     /// @note the weights are calculated as \f$ w = w_x * w_y \f$
489 04 Jan 06 peter 122     ///
489 04 Jan 06 peter 123     /// @return AveragerWeighted for x
489 04 Jan 06 peter 124     ///
718 26 Dec 06 jari 125     const AveragerWeighted& x_averager(void) const;
489 04 Jan 06 peter 126
489 04 Jan 06 peter 127     ///
489 04 Jan 06 peter 128     /// @note the weights are calculated as \f$ w = w_x * w_y \f$
489 04 Jan 06 peter 129     ///
489 04 Jan 06 peter 130     /// @return AveragerWeighted for y
489 04 Jan 06 peter 131     ///
718 26 Dec 06 jari 132     const AveragerWeighted& y_averager(void) const;
489 04 Jan 06 peter 133
2570 26 Sep 11 peter 134     /**
2570 26 Sep 11 peter 135        \brief Addition assignment operator.
2570 26 Sep 11 peter 136
2570 26 Sep 11 peter 137        \return resulting AveragerPairWeighted
2570 26 Sep 11 peter 138
2570 26 Sep 11 peter 139        \since New in yat 0.8
2570 26 Sep 11 peter 140      */
2570 26 Sep 11 peter 141     const AveragerPairWeighted& operator+=(const AveragerPairWeighted&);
2570 26 Sep 11 peter 142
474 22 Dec 05 peter 143   private:
489 04 Jan 06 peter 144     AveragerWeighted x_; // weighted averager with w = w_x*w_y
489 04 Jan 06 peter 145     AveragerWeighted y_; // weighted averager with w = w_x*w_y
2568 26 Sep 11 peter 146     double wxy_centered_;
2568 26 Sep 11 peter 147
2568 26 Sep 11 peter 148     void xy_add(double mx, double my, double xy_centered, double w);
474 22 Dec 05 peter 149   };
474 22 Dec 05 peter 150
890 25 Sep 07 markus 151   /**
915 29 Sep 07 peter 152      \brief adding a ranges of values to AveragerPairWeighted \a ap
1886 31 Mar 09 peter 153
3530 11 Oct 16 peter 154      Type Requirement:
3530 11 Oct 16 peter 155      - \c InputIterator1 models \ref concept_data_iterator
3530 11 Oct 16 peter 156      - \c InputIterator1 is \single_pass_iterator
3530 11 Oct 16 peter 157      - \c InputIterator2 models \ref concept_data_iterator
3530 11 Oct 16 peter 158      - \c InputIterator2 is \single_pass_iterator
3530 11 Oct 16 peter 159
1886 31 Mar 09 peter 160      \relates AveragerPairWeighted
914 29 Sep 07 peter 161    */
2202 21 Feb 10 peter 162   template <class InputIterator1, class InputIterator2>
4200 19 Aug 22 peter 163   void add(AveragerPairWeighted& ap, InputIterator1 first1,
2202 21 Feb 10 peter 164            InputIterator1 last1, InputIterator2 first2)
890 25 Sep 07 markus 165   {
3530 11 Oct 16 peter 166     BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<InputIterator1>));
3530 11 Oct 16 peter 167     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator1>));
3530 11 Oct 16 peter 168     BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<InputIterator2>));
3530 11 Oct 16 peter 169     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator2>));
3530 11 Oct 16 peter 170
2202 21 Feb 10 peter 171     utility::iterator_traits<InputIterator1> traits1;
2202 21 Feb 10 peter 172     utility::iterator_traits<InputIterator2> traits2;
916 30 Sep 07 peter 173     for ( ; first1 != last1; ++first1, ++first2) {
4200 19 Aug 22 peter 174       ap.add(traits1.data(first1), traits2.data(first2),
2202 21 Feb 10 peter 175              traits1.weight(first1), traits2.weight(first2));
916 30 Sep 07 peter 176     }
890 25 Sep 07 markus 177   }
890 25 Sep 07 markus 178
890 25 Sep 07 markus 179
1043 06 Feb 08 peter 180   /**
1043 06 Feb 08 peter 181      \brief adding four ranges of values to AveragerPairWeighted \a ap
1886 31 Mar 09 peter 182
2202 21 Feb 10 peter 183      Iterators must be unweighted.
2202 21 Feb 10 peter 184
3530 11 Oct 16 peter 185      Type Requirement:
3530 11 Oct 16 peter 186      - \c InputIterator1 models \single_pass_iterator
3530 11 Oct 16 peter 187      - \c InputIterator1 models \readable_iterator
3530 11 Oct 16 peter 188      - \c InputIterator1 is weighted
3530 11 Oct 16 peter 189      - \c InputIterator2 models \single_pass_iterator
3530 11 Oct 16 peter 190      - \c InputIterator2 models \readable_iterator
3530 11 Oct 16 peter 191      - \c InputIterator2 is weighted
3530 11 Oct 16 peter 192      - \c InputIterator3 models \single_pass_iterator
3530 11 Oct 16 peter 193      - \c InputIterator3 models \readable_iterator
3530 11 Oct 16 peter 194      - \c InputIterator3 is weighted
3530 11 Oct 16 peter 195      - \c InputIterator4 models \single_pass_iterator
3530 11 Oct 16 peter 196      - \c InputIterator4 models \readable_iterator
3530 11 Oct 16 peter 197      - \c InputIterator4 is weighted
3530 11 Oct 16 peter 198
1886 31 Mar 09 peter 199      \relates AveragerPairWeighted
1043 06 Feb 08 peter 200    */
2202 21 Feb 10 peter 201   template <typename InputIterator1, typename InputIterator2
2202 21 Feb 10 peter 202             , typename InputIterator3, typename InputIterator4>
3530 11 Oct 16 peter 203   void add(AveragerPairWeighted& ap, InputIterator1 x, InputIterator1 xlast,
2202 21 Feb 10 peter 204            InputIterator2 y, InputIterator3 wx, InputIterator4 wy)
582 09 May 06 markus 205   {
3530 11 Oct 16 peter 206     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator1>));
3530 11 Oct 16 peter 207     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator1>));
3530 11 Oct 16 peter 208     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator2>));
3530 11 Oct 16 peter 209     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator2>));
3530 11 Oct 16 peter 210     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator3>));
3530 11 Oct 16 peter 211     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator3>));
3530 11 Oct 16 peter 212     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator4>));
3530 11 Oct 16 peter 213     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator4>));
3530 11 Oct 16 peter 214
1043 06 Feb 08 peter 215     utility::check_iterator_is_unweighted(x);
1043 06 Feb 08 peter 216     utility::check_iterator_is_unweighted(y);
1043 06 Feb 08 peter 217     utility::check_iterator_is_unweighted(wx);
1043 06 Feb 08 peter 218     utility::check_iterator_is_unweighted(wy);
1043 06 Feb 08 peter 219     while (x!=xlast) {
1043 06 Feb 08 peter 220       ap.add(*x, *y, *wx, *wy);
1043 06 Feb 08 peter 221       ++x;
1043 06 Feb 08 peter 222       ++y;
1043 06 Feb 08 peter 223       ++wx;
1043 06 Feb 08 peter 224       ++wy;
916 30 Sep 07 peter 225     }
582 09 May 06 markus 226   }
582 09 May 06 markus 227
1043 06 Feb 08 peter 228
683 11 Oct 06 jari 229 }}} // of namespace statistics, yat, and theplu
582 09 May 06 markus 230
474 22 Dec 05 peter 231 #endif