yat/statistics/AveragerWeighted.h

Code
Comments
Other
Rev Date Author Line
683 11 Oct 06 jari 1 #ifndef _theplu_yat_statistics_averagerweighted_
683 11 Oct 06 jari 2 #define _theplu_yat_statistics_averagerweighted_
94 09 Jun 04 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) 2004 Jari Häkkinen, Peter Johansson, Cecilia Ritz
2119 12 Dec 09 peter 8   Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 9   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
3550 03 Jan 17 peter 10   Copyright (C) 2009, 2010, 2011, 2016 Peter Johansson
94 09 Jun 04 peter 11
1437 25 Aug 08 peter 12   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 13
675 10 Oct 06 jari 14   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 15   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 16   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 17   License, or (at your option) any later version.
675 10 Oct 06 jari 18
675 10 Oct 06 jari 19   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 20   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 21   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 22   General Public License for more details.
675 10 Oct 06 jari 23
675 10 Oct 06 jari 24   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 25   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 26 */
675 10 Oct 06 jari 27
3528 11 Oct 16 peter 28 #include "yat/utility/concept_check.h"
937 05 Oct 07 peter 29 #include "yat/utility/iterator_traits.h"
675 10 Oct 06 jari 30
2202 21 Feb 10 peter 31 #include <boost/concept_check.hpp>
2202 21 Feb 10 peter 32
397 13 Oct 05 jari 33 #include <cmath>
295 29 Apr 05 peter 34
94 09 Jun 04 peter 35 namespace theplu{
680 11 Oct 06 jari 36 namespace yat{
197 27 Oct 04 jari 37 namespace statistics{
295 29 Apr 05 peter 38
94 09 Jun 04 peter 39   ///
486 04 Jan 06 peter 40   /// @brief Class to calulate averages with weights.
94 09 Jun 04 peter 41   ///
486 04 Jan 06 peter 42   /// There are several different reasons why a statistical analysis
486 04 Jan 06 peter 43   /// needs to adjust for weighting. In the litterature reasons are
486 04 Jan 06 peter 44   /// mainly divided into two kinds of weights - probablity weights
486 04 Jan 06 peter 45   /// and analytical weights. 1) Analytical weights are appropriate
486 04 Jan 06 peter 46   /// for scientific experiments where some measurements are known to
486 04 Jan 06 peter 47   /// be more precise than others. The larger weight a measurement has
486 04 Jan 06 peter 48   /// the more precise is is assumed to be, or more formally the
3528 11 Oct 16 peter 49   /// weight is proportional to the reciprocal variance
589 24 Aug 06 peter 50   /// \f$ \sigma_i^2 = \frac{\sigma^2}{w_i} \f$. 2) Probability weights
486 04 Jan 06 peter 51   /// are used for the situation when calculating averages over a
888 24 Sep 07 peter 52   /// distribution \f$ f \f$ , but sampling from a distribution \f$ f' \f$.
888 24 Sep 07 peter 53   /// Compensating for this discrepancy averages of observables
486 04 Jan 06 peter 54   /// are taken to be \f$ \sum \frac{f}{f'}X \f$ For further discussion:
1185 28 Feb 08 peter 55   /// see \ref weighted_statistics
3528 11 Oct 16 peter 56   ///
486 04 Jan 06 peter 57   /// If nothing else stated, each function fulfills the
486 04 Jan 06 peter 58   /// following:<br> <ul><li>Setting a weight to zero corresponds to
486 04 Jan 06 peter 59   /// removing the data point from the dataset.</li><li> Setting all
486 04 Jan 06 peter 60   /// weights to unity, the yields the same result as from
486 04 Jan 06 peter 61   /// corresponding function in Averager.</li><li> Rescaling weights
486 04 Jan 06 peter 62   /// does not change the performance of the object.</li></ul>
194 27 Oct 04 jari 63   ///
486 04 Jan 06 peter 64   /// @see Averager AveragerPair AveragerPairWeighted
486 04 Jan 06 peter 65   ///
295 29 Apr 05 peter 66   class AveragerWeighted
94 09 Jun 04 peter 67   {
94 09 Jun 04 peter 68   public:
94 09 Jun 04 peter 69
94 09 Jun 04 peter 70     ///
703 18 Dec 06 jari 71     /// @brief The default constructor
3528 11 Oct 16 peter 72     ///
703 18 Dec 06 jari 73     AveragerWeighted(void);
95 09 Jun 04 peter 74
94 09 Jun 04 peter 75     ///
703 18 Dec 06 jari 76     /// Adding a data point \a d, with weight \a w (default is 1)
94 09 Jun 04 peter 77     ///
703 18 Dec 06 jari 78     void  add(const double d,const double w=1);
94 09 Jun 04 peter 79
94 09 Jun 04 peter 80     ///
3528 11 Oct 16 peter 81     /// @brief Calculate the weighted mean
94 09 Jun 04 peter 82     ///
3528 11 Oct 16 peter 83     /// @return \f$ \frac{\sum w_ix_i}{\sum w_i} \f$
94 09 Jun 04 peter 84     ///
718 26 Dec 06 jari 85     double mean(void) const;
699 26 Oct 06 peter 86
94 09 Jun 04 peter 87     ///
3528 11 Oct 16 peter 88     /// @brief Weighted version of number of data points.
699 26 Oct 06 peter 89     ///
699 26 Oct 06 peter 90     /// If all
589 24 Aug 06 peter 91     /// weights are equal, the unweighted version is identical to the
589 24 Aug 06 peter 92     /// non-weighted version. Adding a data point with zero weight
589 24 Aug 06 peter 93     /// does not change n(). The calculated value is always smaller
589 24 Aug 06 peter 94     /// than the actual number of data points added to object.
589 24 Aug 06 peter 95     ///
492 09 Jan 06 peter 96     /// @return \f$ \frac{\left(\sum w_i\right)^2}{\sum w_i^2} \f$
3528 11 Oct 16 peter 97     ///
718 26 Dec 06 jari 98     double n(void) const;
492 09 Jan 06 peter 99
492 09 Jan 06 peter 100     ///
703 18 Dec 06 jari 101     /// @brief Rescale object.
197 27 Oct 04 jari 102     ///
703 18 Dec 06 jari 103     /// Each data point is rescaled as \f$ x = a * x \f$
703 18 Dec 06 jari 104     ///
703 18 Dec 06 jari 105     void rescale(double a);
197 27 Oct 04 jari 106
197 27 Oct 04 jari 107     ///
703 18 Dec 06 jari 108     /// @brief Reset everything to zero.
96 09 Jun 04 peter 109     ///
703 18 Dec 06 jari 110     void reset(void);
96 09 Jun 04 peter 111
96 09 Jun 04 peter 112     ///
703 18 Dec 06 jari 113     /// @brief The standard deviation is defined as the square root of
703 18 Dec 06 jari 114     /// the variance().
118 20 Jul 04 peter 115     ///
486 04 Jan 06 peter 116     /// @return The standard deviation, root of the variance().
200 01 Nov 04 peter 117     ///
718 26 Dec 06 jari 118     double std(void) const;
200 01 Nov 04 peter 119
200 01 Nov 04 peter 120     ///
703 18 Dec 06 jari 121     /// @brief Calculates standard deviation of the mean().
703 18 Dec 06 jari 122     ///
703 18 Dec 06 jari 123     /// Variance from the
486 04 Jan 06 peter 124     /// weights are here neglected. This is true when the weight is
486 04 Jan 06 peter 125     /// known before the measurement. In case this is not a good
486 04 Jan 06 peter 126     /// approximation, use bootstrapping to estimate the error.
397 13 Oct 05 jari 127     ///
3528 11 Oct 16 peter 128     /// @return \f$ \frac{\sum w^2}{\left(\sum w\right)^3}\sum w(x-m)^2 \f$
3528 11 Oct 16 peter 129     /// where \f$ m \f$ is the mean()
397 13 Oct 05 jari 130     ///
718 26 Dec 06 jari 131     double standard_error(void) const;
397 13 Oct 05 jari 132
397 13 Oct 05 jari 133     ///
3528 11 Oct 16 peter 134     /// Calculating the sum of weights
95 09 Jun 04 peter 135     ///
486 04 Jan 06 peter 136     /// @return \f$ \sum w_i \f$
486 04 Jan 06 peter 137     ///
718 26 Dec 06 jari 138     double sum_w(void) const;
95 09 Jun 04 peter 139
95 09 Jun 04 peter 140     ///
492 09 Jan 06 peter 141     /// @return \f$ \sum w_i^2 \f$
492 09 Jan 06 peter 142     ///
718 26 Dec 06 jari 143     double sum_ww(void) const;
492 09 Jan 06 peter 144
492 09 Jan 06 peter 145     ///
3528 11 Oct 16 peter 146     /// \f$ \sum w_ix_i \f$
489 04 Jan 06 peter 147     ///
489 04 Jan 06 peter 148     /// @return weighted sum of x
489 04 Jan 06 peter 149     ///
718 26 Dec 06 jari 150     double sum_wx(void) const;
489 04 Jan 06 peter 151
489 04 Jan 06 peter 152     ///
772 25 Feb 07 peter 153     ///  @return \f$ \sum w_i x_i^2 \f$
772 25 Feb 07 peter 154     ///
772 25 Feb 07 peter 155     double sum_wxx(void) const;
772 25 Feb 07 peter 156
772 25 Feb 07 peter 157     ///
486 04 Jan 06 peter 158     /// @return \f$ \sum_i w_i (x_i-m)^2\f$
475 22 Dec 05 peter 159     ///
718 26 Dec 06 jari 160     double sum_xx_centered(void) const;
475 22 Dec 05 peter 161
699 26 Oct 06 peter 162     /**
699 26 Oct 06 peter 163        The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
3528 11 Oct 16 peter 164        }{\sum w_i} \f$, where \a m is the known mean.
3528 11 Oct 16 peter 165
699 26 Oct 06 peter 166        @return Variance when the mean is known to be \a m.
699 26 Oct 06 peter 167     */
718 26 Dec 06 jari 168     double variance(const double m) const;
200 01 Nov 04 peter 169
699 26 Oct 06 peter 170     /**
699 26 Oct 06 peter 171        The variance is calculated as \f$ \frac{\sum w_i (x_i - m)^2
699 26 Oct 06 peter 172        }{\sum w_i} \f$, where \a m is the mean(). Here the weight are
699 26 Oct 06 peter 173        interpreted as probability weights. For analytical weights the
699 26 Oct 06 peter 174        variance has no meaning as each data point has its own
699 26 Oct 06 peter 175        variance.
3528 11 Oct 16 peter 176
699 26 Oct 06 peter 177        @return The variance.
699 26 Oct 06 peter 178     */
718 26 Dec 06 jari 179     double variance(void) const;
400 18 Oct 05 peter 180
2032 19 Aug 09 peter 181     /**
2032 19 Aug 09 peter 182        operator to add an AveragerWeighted
414 01 Dec 05 peter 183
2032 19 Aug 09 peter 184        \since New in yat 0.6
2032 19 Aug 09 peter 185     */
2032 19 Aug 09 peter 186     const AveragerWeighted& operator+=(const AveragerWeighted&);
3528 11 Oct 16 peter 187
400 18 Oct 05 peter 188   private:
95 09 Jun 04 peter 189     ///
718 26 Dec 06 jari 190     ///  @return \f$ \sum w_i^2x_i \f$
718 26 Dec 06 jari 191     ///
718 26 Dec 06 jari 192     double sum_wwx(void) const;
718 26 Dec 06 jari 193
718 26 Dec 06 jari 194     ///
95 09 Jun 04 peter 195     ///  @return \f$ \sum w_i^2x_i^2 \f$
95 09 Jun 04 peter 196     ///
718 26 Dec 06 jari 197     double sum_wwxx(void) const;
3528 11 Oct 16 peter 198
2562 25 Sep 11 peter 199     double mean_; // wx/w
2562 25 Sep 11 peter 200     double m2_; // wxx - m*m*w
2562 25 Sep 11 peter 201     double w_;
2562 25 Sep 11 peter 202     double ww_;
2562 25 Sep 11 peter 203     double wwxx_;
2562 25 Sep 11 peter 204     double wxx_;
718 26 Dec 06 jari 205
2562 25 Sep 11 peter 206     // double wwx; // m2_ + m*m*w
2562 25 Sep 11 peter 207     // double wx; // w*mean_
94 09 Jun 04 peter 208   };
94 09 Jun 04 peter 209
910 29 Sep 07 peter 210   /**
916 30 Sep 07 peter 211      \brief adding a range of values to AveragerWeighted \a a
1043 06 Feb 08 peter 212
2202 21 Feb 10 peter 213      If InputIterator is non-weighted unitary weights are used.
1886 31 Mar 09 peter 214
3528 11 Oct 16 peter 215      Type Requirement:
3528 11 Oct 16 peter 216      - \c InputIterator models \ref concept_data_iterator
3528 11 Oct 16 peter 217      - \c InputIterator is \single_pass_iterator
3528 11 Oct 16 peter 218
1886 31 Mar 09 peter 219      \relates AveragerWeighted
910 29 Sep 07 peter 220    */
2202 21 Feb 10 peter 221   template <typename InputIterator>
2202 21 Feb 10 peter 222   void add(AveragerWeighted& a, InputIterator first, InputIterator last)
910 29 Sep 07 peter 223   {
3528 11 Oct 16 peter 224     BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<InputIterator>));
3528 11 Oct 16 peter 225     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator>));
2202 21 Feb 10 peter 226     utility::iterator_traits<InputIterator> traits;
910 29 Sep 07 peter 227     for ( ; first != last; ++first)
2159 19 Jan 10 peter 228       a.add(traits.data(first), traits.weight(first));
910 29 Sep 07 peter 229   }
910 29 Sep 07 peter 230
916 30 Sep 07 peter 231   /**
3528 11 Oct 16 peter 232      \brief add values from two ranges to AveragerWeighted \a a
916 30 Sep 07 peter 233
916 30 Sep 07 peter 234      Add data from range [first1, last1) with their corresponding
916 30 Sep 07 peter 235      weight in range [first2, first2 + distance(first, last) ).
916 30 Sep 07 peter 236
3528 11 Oct 16 peter 237      Type Requirement:
3528 11 Oct 16 peter 238      - \c InputIterator1 models \single_pass_iterator
3528 11 Oct 16 peter 239      - \c InputIterator1 models \readable_iterator
3528 11 Oct 16 peter 240      - \c InputIterator1 is unweighted
3528 11 Oct 16 peter 241      - \c InputIterator2 models \single_pass_iterator
3528 11 Oct 16 peter 242      - \c InputIterator2 models \readable_iterator
3528 11 Oct 16 peter 243      - \c InputIterator2 is unweighted
1886 31 Mar 09 peter 244
1886 31 Mar 09 peter 245      \relates AveragerWeighted
916 30 Sep 07 peter 246    */
2202 21 Feb 10 peter 247   template <typename InputIterator1, typename InputIterator2>
3528 11 Oct 16 peter 248   void add(AveragerWeighted& a, InputIterator1 first1, InputIterator1 last1,
2202 21 Feb 10 peter 249            InputIterator2 first2)
582 09 May 06 markus 250   {
3528 11 Oct 16 peter 251     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator1>));
3528 11 Oct 16 peter 252     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator1>));
3528 11 Oct 16 peter 253     BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<InputIterator2>));
3528 11 Oct 16 peter 254     BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<InputIterator2>));
3528 11 Oct 16 peter 255
916 30 Sep 07 peter 256     utility::check_iterator_is_unweighted(first1);
916 30 Sep 07 peter 257     utility::check_iterator_is_unweighted(first2);
916 30 Sep 07 peter 258     for ( ; first1 != last1; ++first1, ++first2)
916 30 Sep 07 peter 259       a.add(*first1, *first2);
582 09 May 06 markus 260   }
582 09 May 06 markus 261
683 11 Oct 06 jari 262 }}} // of namespace statistics, yat, and theplu
94 09 Jun 04 peter 263
94 09 Jun 04 peter 264 #endif

Generated on Sat Oct 14 02:23:39 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)