yat/statistics/averager_base.h

Code
Comments
Other
Rev Date Author Line
2809 06 Aug 12 peter 1 #ifndef _theplu_yat_statistics_average_base
2809 06 Aug 12 peter 2 #define _theplu_yat_statistics_average_base
2809 06 Aug 12 peter 3
2809 06 Aug 12 peter 4 // $Id$
2809 06 Aug 12 peter 5
2809 06 Aug 12 peter 6 /*
2809 06 Aug 12 peter 7   Copyright (C) 2012 Peter Johansson
2809 06 Aug 12 peter 8
2809 06 Aug 12 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
2809 06 Aug 12 peter 10
2809 06 Aug 12 peter 11   The yat library is free software; you can redistribute it and/or
2809 06 Aug 12 peter 12   modify it under the terms of the GNU General Public License as
2809 06 Aug 12 peter 13   published by the Free Software Foundation; either version 3 of the
2809 06 Aug 12 peter 14   License, or (at your option) any later version.
2809 06 Aug 12 peter 15
2809 06 Aug 12 peter 16   The yat library is distributed in the hope that it will be useful,
2809 06 Aug 12 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
2809 06 Aug 12 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2809 06 Aug 12 peter 19   General Public License for more details.
2809 06 Aug 12 peter 20
2809 06 Aug 12 peter 21   You should have received a copy of the GNU General Public License
2809 06 Aug 12 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
2809 06 Aug 12 peter 23 */
2809 06 Aug 12 peter 24
2809 06 Aug 12 peter 25 #include <boost/concept_check.hpp>
2809 06 Aug 12 peter 26
2809 06 Aug 12 peter 27 #include <cmath>
2809 06 Aug 12 peter 28 #include <limits>
2809 06 Aug 12 peter 29
2809 06 Aug 12 peter 30 namespace theplu {
2809 06 Aug 12 peter 31 namespace yat {
2809 06 Aug 12 peter 32 namespace statistics {
2809 06 Aug 12 peter 33
2809 06 Aug 12 peter 34   /**
2809 06 Aug 12 peter 35      \brief Base class for averager classes
2809 06 Aug 12 peter 36
2809 06 Aug 12 peter 37      This is a base class templated on the derived class (see
2809 06 Aug 12 peter 38      curiously recurring template pattern), which implies there is no
2809 06 Aug 12 peter 39      dynamic polymorphism but the concrete type must be decided at
2809 06 Aug 12 peter 40      compile time.
2809 06 Aug 12 peter 41
2809 06 Aug 12 peter 42      The derived class, \c MyClass, needs to have the following features:
2809 06 Aug 12 peter 43      - \c <a href=http://www.sgi.com/tech/stl/DefaultConstructible.html>
2809 06 Aug 12 peter 44           DefaultConstructible</a>
2809 06 Aug 12 peter 45      - \c <a href=\boost_url/utility/Assignable.html> Assignable</a>
2809 06 Aug 12 peter 46      - \c a function add_impl(double, long int)
2809 06 Aug 12 peter 47      - \c a function rescale_impl(double)
2809 06 Aug 12 peter 48
2809 06 Aug 12 peter 49      These two functions need to be public or if you dont want to
2809 06 Aug 12 peter 50      expose them to the world, make them private and declare \c
2809 06 Aug 12 peter 51      averager_base<Derived> as friend.
2809 06 Aug 12 peter 52
2809 06 Aug 12 peter 53      A derived class is typically declared as
2809 06 Aug 12 peter 54
2809 06 Aug 12 peter 55      \code
2809 06 Aug 12 peter 56      class MyClass : public averager_base<MyClass>
2809 06 Aug 12 peter 57      {
2809 06 Aug 12 peter 58      ...
2809 06 Aug 12 peter 59      private:
2809 06 Aug 12 peter 60        friend class averager_base<MyClass>;
2809 06 Aug 12 peter 61        void add_impl(double, long int);
2809 06 Aug 12 peter 62        void rescale_impl(double);
2809 06 Aug 12 peter 63      };
2809 06 Aug 12 peter 64      \endcode
2809 06 Aug 12 peter 65
2809 06 Aug 12 peter 66      \since yat 0.9
2809 06 Aug 12 peter 67    */
2809 06 Aug 12 peter 68   template<class Derived>
2809 06 Aug 12 peter 69   class averager_base
2809 06 Aug 12 peter 70   {
2809 06 Aug 12 peter 71   public:
2809 06 Aug 12 peter 72     /**
2809 06 Aug 12 peter 73        \brief Constructor
2809 06 Aug 12 peter 74      */
2809 06 Aug 12 peter 75     averager_base(void) : n_(0), mean_(0)
2809 06 Aug 12 peter 76     {
2809 06 Aug 12 peter 77       BOOST_CONCEPT_ASSERT((boost::DefaultConstructible<Derived>));
2809 06 Aug 12 peter 78     }
2809 06 Aug 12 peter 79
2809 06 Aug 12 peter 80     /**
2809 06 Aug 12 peter 81        \brief Destructor
2809 06 Aug 12 peter 82      */
2809 06 Aug 12 peter 83     virtual ~averager_base(void)=0;
2809 06 Aug 12 peter 84
2809 06 Aug 12 peter 85     /**
2809 06 Aug 12 peter 86        \brief add a data point
2809 06 Aug 12 peter 87
2809 06 Aug 12 peter 88        Adding \a n number of data point(s) with value \a x.
2809 06 Aug 12 peter 89     */
2809 06 Aug 12 peter 90     void add(double x, long n=1) {static_cast<Derived*>(this)->add_impl(x, n);}
2809 06 Aug 12 peter 91
2809 06 Aug 12 peter 92     /**
2809 06 Aug 12 peter 93        \brief mean
2809 06 Aug 12 peter 94
2809 06 Aug 12 peter 95        \return %Mean of presented data, \f$ \frac{1}{n}\sum x_i \f$
2809 06 Aug 12 peter 96     */
2809 06 Aug 12 peter 97     double mean(void) const
2809 06 Aug 12 peter 98     { return n_ ? mean_ : std::numeric_limits<double>::quiet_NaN(); }
2809 06 Aug 12 peter 99
2809 06 Aug 12 peter 100     /**
2809 06 Aug 12 peter 101        \brief number of data points
2809 06 Aug 12 peter 102
2809 06 Aug 12 peter 103        \return Number of data points
2809 06 Aug 12 peter 104     */
2809 06 Aug 12 peter 105     long n(void) const { return n_; }
2809 06 Aug 12 peter 106
2809 06 Aug 12 peter 107     /**
2809 06 Aug 12 peter 108        \brief Rescales the object
2809 06 Aug 12 peter 109
2809 06 Aug 12 peter 110        \f$ \forall x_i \rightarrow a*x_i \f$,
2809 06 Aug 12 peter 111     */
2809 06 Aug 12 peter 112     void rescale(double a) { static_cast<Derived*>(this)->rescale_impl(a); }
2809 06 Aug 12 peter 113
2809 06 Aug 12 peter 114     /**
2809 06 Aug 12 peter 115        \brief Reset object
2809 06 Aug 12 peter 116
2809 06 Aug 12 peter 117        Restore Averager as if data were never added
2809 06 Aug 12 peter 118     */
2809 06 Aug 12 peter 119     void reset(void)
2809 06 Aug 12 peter 120     {
2809 06 Aug 12 peter 121       Derived tmp;
2809 06 Aug 12 peter 122       *static_cast<Derived*>(this) = tmp;
2809 06 Aug 12 peter 123     }
2809 06 Aug 12 peter 124
2809 06 Aug 12 peter 125     /**
2809 06 Aug 12 peter 126        \return The sum of data values
2809 06 Aug 12 peter 127     */
2809 06 Aug 12 peter 128     double sum_x(void)  const { return n_*mean_; }
2809 06 Aug 12 peter 129
2809 06 Aug 12 peter 130   protected:
2809 06 Aug 12 peter 131     /**
2809 06 Aug 12 peter 132        number of data points
2809 06 Aug 12 peter 133     */
2809 06 Aug 12 peter 134     long int n_;
2809 06 Aug 12 peter 135
2809 06 Aug 12 peter 136     /**
2809 06 Aug 12 peter 137        \brief mean
2809 06 Aug 12 peter 138     */
2809 06 Aug 12 peter 139     double mean_;
2809 06 Aug 12 peter 140
2809 06 Aug 12 peter 141     /**
2809 06 Aug 12 peter 142        add \a n data points with value \a x
2809 06 Aug 12 peter 143     */
2809 06 Aug 12 peter 144     void add1(double x, long int n);
2809 06 Aug 12 peter 145     /**
2809 06 Aug 12 peter 146        add one data point with value \a delta + mean()
2809 06 Aug 12 peter 147      */
2809 06 Aug 12 peter 148     void add1(double delta);
2809 06 Aug 12 peter 149
2809 06 Aug 12 peter 150     /**
2809 06 Aug 12 peter 151        rescale as \f$ mean \rightarrow x mean \f$
2809 06 Aug 12 peter 152      */
2809 06 Aug 12 peter 153     double rescale1(double x);
2809 06 Aug 12 peter 154   };
2809 06 Aug 12 peter 155
2809 06 Aug 12 peter 156
2809 06 Aug 12 peter 157   /**
2809 06 Aug 12 peter 158      \brief Base class for averagers calculating mean and variance.
2809 06 Aug 12 peter 159
2809 06 Aug 12 peter 160      This is an extension of averager_base and provides more
2809 06 Aug 12 peter 161      functionality related to the second moment (variance).
2809 06 Aug 12 peter 162
2809 06 Aug 12 peter 163      For requirement on \c Derived see averager_base.
2809 06 Aug 12 peter 164
2809 06 Aug 12 peter 165      \see Averager which is an example derived class
2809 06 Aug 12 peter 166
2809 06 Aug 12 peter 167      \since yat 0.9
2809 06 Aug 12 peter 168    */
2809 06 Aug 12 peter 169   template<class Derived>
2809 06 Aug 12 peter 170   class averager_base2 : public averager_base<Derived>
2809 06 Aug 12 peter 171   {
2809 06 Aug 12 peter 172   public:
2809 06 Aug 12 peter 173     /**
2809 06 Aug 12 peter 174        \brief Constructor
2809 06 Aug 12 peter 175      */
2809 06 Aug 12 peter 176     averager_base2(void) : cm2_(0) {}
2809 06 Aug 12 peter 177
2809 06 Aug 12 peter 178     /**
2809 06 Aug 12 peter 179        \brief Destructor
2809 06 Aug 12 peter 180     */
2809 06 Aug 12 peter 181     virtual ~averager_base2(void)=0;
2809 06 Aug 12 peter 182
2809 06 Aug 12 peter 183     /**
2809 06 Aug 12 peter 184        \brief Coeffient of variation
2809 06 Aug 12 peter 185
2809 06 Aug 12 peter 186        Coeffient of variation (cv) is defined as ratio between the
2809 06 Aug 12 peter 187        standard deviation and the mean: \f$ \frac{\sigma}{\mu} \f$.
2809 06 Aug 12 peter 188
2809 06 Aug 12 peter 189        @return standard deviation divided by mean.
2809 06 Aug 12 peter 190     */
2809 06 Aug 12 peter 191     double cv(void) const { return std() / this->mean(); }
2809 06 Aug 12 peter 192
2809 06 Aug 12 peter 193     /**
2809 06 Aug 12 peter 194        \return Standard error, i.e. standard deviation of the mean
2809 06 Aug 12 peter 195        \f$ \sqrt{variance()/n} \f$
2809 06 Aug 12 peter 196     */
2809 06 Aug 12 peter 197     double standard_error(void) const
2809 06 Aug 12 peter 198     { return std::sqrt(variance()/this->n()); }
2809 06 Aug 12 peter 199
2809 06 Aug 12 peter 200     /**
2809 06 Aug 12 peter 201        \brief The standard deviation is defined as the square root of
2809 06 Aug 12 peter 202        the variance.
2809 06 Aug 12 peter 203
2809 06 Aug 12 peter 204        \return The standard deviation, root of the variance().
2809 06 Aug 12 peter 205     */
2809 06 Aug 12 peter 206     double std(void) const { return std::sqrt(variance()); }
2809 06 Aug 12 peter 207
2809 06 Aug 12 peter 208     /**
2809 06 Aug 12 peter 209        \brief The standard deviation is defined as the square root of
2809 06 Aug 12 peter 210        the variance.
2809 06 Aug 12 peter 211
2809 06 Aug 12 peter 212        \return Standard deviation around \a m, root of the variance(m).
2809 06 Aug 12 peter 213     */
2809 06 Aug 12 peter 214     double std(double m) const { return std::sqrt(variance(m)); }
2809 06 Aug 12 peter 215
2809 06 Aug 12 peter 216     /**
2809 06 Aug 12 peter 217        \return The sum of squares
2809 06 Aug 12 peter 218     */
2809 06 Aug 12 peter 219     double sum_xx(void) const
2809 06 Aug 12 peter 220     { return cm2_ + this->n_*this->mean_*this->mean_; }
2809 06 Aug 12 peter 221
2809 06 Aug 12 peter 222     /**
2809 06 Aug 12 peter 223        \return \f$ \sum_i (x_i-m)^2 \f$
2809 06 Aug 12 peter 224     */
2809 06 Aug 12 peter 225     double sum_xx_centered(void)  const { return cm2_; }
2809 06 Aug 12 peter 226
2809 06 Aug 12 peter 227     /**
2809 06 Aug 12 peter 228        @brief The variance with known mean
2809 06 Aug 12 peter 229
2809 06 Aug 12 peter 230        The variance is calculated as
2809 06 Aug 12 peter 231        \f$ \frac{1}{n}\sum (x_i-m)^2 \f$.
2809 06 Aug 12 peter 232
2809 06 Aug 12 peter 233        \return Variance when the mean is known to be \a m.
2809 06 Aug 12 peter 234     */
2809 06 Aug 12 peter 235     double variance(double m) const
2809 06 Aug 12 peter 236     { return sum_xx()/this->n() + m*m - 2*m*this->mean();}
2809 06 Aug 12 peter 237
2809 06 Aug 12 peter 238     /**
2809 06 Aug 12 peter 239        \brief The estimated variance
2809 06 Aug 12 peter 240
2809 06 Aug 12 peter 241        The variance is calculated as \f$ \frac{1}{N}\sum_i
2809 06 Aug 12 peter 242        (x_i-m)^2 \f$, where \f$ m \f$ is the mean.
2809 06 Aug 12 peter 243
2809 06 Aug 12 peter 244        \return Estimation of variance
2809 06 Aug 12 peter 245     */
2809 06 Aug 12 peter 246     double variance(void) const { return sum_xx_centered()/this->n_; }
2809 06 Aug 12 peter 247
2809 06 Aug 12 peter 248     /**
2809 06 Aug 12 peter 249        The variance is calculated using the \f$ (n-1) \f$ correction,
2809 06 Aug 12 peter 250        which means it is the best unbiased estimator of the variance
2809 06 Aug 12 peter 251        \f$ \frac{1}{N-1}\sum_i (x_i-m)^2 \f$, where \f$ m \f$ is the
2809 06 Aug 12 peter 252        mean.
2809 06 Aug 12 peter 253
2809 06 Aug 12 peter 254        @return unbiased estimation of variance
2809 06 Aug 12 peter 255     */
2809 06 Aug 12 peter 256     double variance_unbiased(void) const
2809 06 Aug 12 peter 257     {
2809 06 Aug 12 peter 258       return (this->n_>1) ? sum_xx_centered()/(this->n_-1) :
2809 06 Aug 12 peter 259         std::numeric_limits<double>::quiet_NaN();
2809 06 Aug 12 peter 260     }
2809 06 Aug 12 peter 261
2809 06 Aug 12 peter 262   protected:
2809 06 Aug 12 peter 263     /**
2809 06 Aug 12 peter 264        sum of values squared values centralized \f$ \sum (x-m)^2 \f$
2809 06 Aug 12 peter 265      */
2809 06 Aug 12 peter 266     double cm2_;
2809 06 Aug 12 peter 267
2809 06 Aug 12 peter 268     /**
2809 06 Aug 12 peter 269        add a set of \a n data points with mean \a mean and centralized
2809 06 Aug 12 peter 270        squaed sum \a cm
2809 06 Aug 12 peter 271      */
2809 06 Aug 12 peter 272     void add2(double mean, double cm2, long int n);
2809 06 Aug 12 peter 273
2809 06 Aug 12 peter 274     /**
2809 06 Aug 12 peter 275        add one data point with value \a delta + mean()
2809 06 Aug 12 peter 276      */
2809 06 Aug 12 peter 277     void add2(double delta);
2809 06 Aug 12 peter 278
2809 06 Aug 12 peter 279     /**
2809 06 Aug 12 peter 280        rescales \f$ cm2 \rightarrow x^2 cm2 \f$ and calles rescale1(double)
2809 06 Aug 12 peter 281      */
2809 06 Aug 12 peter 282     double rescale2(double x);
2809 06 Aug 12 peter 283   private:
2809 06 Aug 12 peter 284   };
2809 06 Aug 12 peter 285
2809 06 Aug 12 peter 286
2809 06 Aug 12 peter 287   /**
2809 06 Aug 12 peter 288      Base class for averagers calculating first (mean), second
2809 06 Aug 12 peter 289      (variance) and third moment.
2809 06 Aug 12 peter 290
2809 06 Aug 12 peter 291      For requirement on \c Derived see averager_base.
2809 06 Aug 12 peter 292
2809 06 Aug 12 peter 293      \since yat 0.9
2809 06 Aug 12 peter 294    */
2809 06 Aug 12 peter 295   template<class Derived>
2809 06 Aug 12 peter 296   class averager_base3 : public averager_base2<Derived>
2809 06 Aug 12 peter 297   {
2809 06 Aug 12 peter 298   public:
2809 06 Aug 12 peter 299     /**
2809 06 Aug 12 peter 300        \brief Constructor
2809 06 Aug 12 peter 301     */
2809 06 Aug 12 peter 302     averager_base3(void) : averager_base2<Derived>(), cm3_(0) {}
2809 06 Aug 12 peter 303
2809 06 Aug 12 peter 304     /**
2809 06 Aug 12 peter 305        \brief Destructor
2809 06 Aug 12 peter 306      */
2809 06 Aug 12 peter 307     virtual ~averager_base3(void)=0;
2809 06 Aug 12 peter 308
2809 06 Aug 12 peter 309     /**
2809 06 Aug 12 peter 310        \brief Third central moment
2809 06 Aug 12 peter 311
2809 06 Aug 12 peter 312        Third central moment is calculated as \f$ \frac{1}{n} \sum
2809 06 Aug 12 peter 313        (x_i-m)^3 \f$ where \c m is the mean.
2809 06 Aug 12 peter 314      */
2809 06 Aug 12 peter 315     double central_moment3(void) const { return cm3_/this->n_; }
2809 06 Aug 12 peter 316
2809 06 Aug 12 peter 317     /**
2809 06 Aug 12 peter 318        \brief Skewness
2809 06 Aug 12 peter 319
2809 06 Aug 12 peter 320        Skewness is calculated as \f$ \frac{1}{n} \sum (x_i-m) /
2809 06 Aug 12 peter 321        \sigma^3 \f$ where \f$ \sigma \f$ is standard deviation.
2809 06 Aug 12 peter 322
2809 06 Aug 12 peter 323        \note this function uses a biased variance estimate [1/n]
2809 06 Aug 12 peter 324        compared to skewness(utility::VectorBase&), which uses an
2809 06 Aug 12 peter 325        unbiased variance estimate [1/(n-1)] and may thus give
2809 06 Aug 12 peter 326        different results for small \c n.
2809 06 Aug 12 peter 327      */
2809 06 Aug 12 peter 328     double skewness(void) const
2809 06 Aug 12 peter 329     { return central_moment3()/std::pow(this->std(),3); }
2809 06 Aug 12 peter 330
2809 06 Aug 12 peter 331     /**
2809 06 Aug 12 peter 332        \return sum cubed centered
2809 06 Aug 12 peter 333
2809 06 Aug 12 peter 334        \f$ \sum (x_i-m)^3 \f$
2809 06 Aug 12 peter 335      */
2809 06 Aug 12 peter 336     double sum_x3_centered(void) const { return cm3_; }
2809 06 Aug 12 peter 337
2809 06 Aug 12 peter 338   protected:
2809 06 Aug 12 peter 339     /**
2809 06 Aug 12 peter 340        sum cubed centered
2809 06 Aug 12 peter 341      */
2809 06 Aug 12 peter 342     double cm3_;
2809 06 Aug 12 peter 343
2809 06 Aug 12 peter 344     /**
2809 06 Aug 12 peter 345        add a set of \a n data points with mean \c mean, sum squared
2809 06 Aug 12 peter 346        centered \a cm2, and sum cubed centered \a cm3.
2809 06 Aug 12 peter 347      */
2809 06 Aug 12 peter 348     void add3(double mean, double cm2, double cm3, long int n);
2809 06 Aug 12 peter 349
2809 06 Aug 12 peter 350     /**
2809 06 Aug 12 peter 351        add one data point with value \c delta + mean()
2809 06 Aug 12 peter 352      */
2809 06 Aug 12 peter 353     void add3(double delta);
2809 06 Aug 12 peter 354
2809 06 Aug 12 peter 355     /**
2809 06 Aug 12 peter 356        rescales \f$ cm3 = x^3 cm3 \f$ and calles rescale2(double)
2809 06 Aug 12 peter 357      */
2809 06 Aug 12 peter 358     double rescale3(double x);
2809 06 Aug 12 peter 359   };
2809 06 Aug 12 peter 360
2809 06 Aug 12 peter 361   /**
2809 06 Aug 12 peter 362      Base class for averagers calculating first (mean), second
2809 06 Aug 12 peter 363      (variance), third moment, and fourth moment.
2809 06 Aug 12 peter 364
2809 06 Aug 12 peter 365      For requirement on \c Derived see averager_base.
2809 06 Aug 12 peter 366
2809 06 Aug 12 peter 367      \since yat 0.9
2809 06 Aug 12 peter 368   */
2809 06 Aug 12 peter 369   template<class Derived>
2809 06 Aug 12 peter 370   class averager_base4 : public averager_base3<Derived>
2809 06 Aug 12 peter 371   {
2809 06 Aug 12 peter 372   public:
2809 06 Aug 12 peter 373     /**
2809 06 Aug 12 peter 374        \brief Constructor
2809 06 Aug 12 peter 375      */
2809 06 Aug 12 peter 376     averager_base4(void) : averager_base3<Derived>(), cm4_(0) {}
2809 06 Aug 12 peter 377
2809 06 Aug 12 peter 378     /**
2809 06 Aug 12 peter 379        \brief Destructor
2809 06 Aug 12 peter 380      */
2809 06 Aug 12 peter 381     virtual ~averager_base4(void)=0;
2809 06 Aug 12 peter 382
2809 06 Aug 12 peter 383     /**
2809 06 Aug 12 peter 384        \brief Fourth central moment
2809 06 Aug 12 peter 385
2809 06 Aug 12 peter 386        Fourth central moment is calculated as \f$ \frac{1}{n} \sum
2809 06 Aug 12 peter 387        (x_i-m)^4 \f$ where \c m is the mean.
2809 06 Aug 12 peter 388      */
2809 06 Aug 12 peter 389     double central_moment4(void) const { return cm4_/this->n_; }
2809 06 Aug 12 peter 390
2809 06 Aug 12 peter 391     /**
2809 06 Aug 12 peter 392        \brief Kurtosis
2809 06 Aug 12 peter 393
2809 06 Aug 12 peter 394        Kurtosis is calculated as \f$ \frac{1}{n} \sum (x-m)^4 /
2809 06 Aug 12 peter 395        \sigma^4 - 3.0 \f$ where \c m is mean, and \f$ \sigma \f$ is
2809 06 Aug 12 peter 396        standard deviation.
2809 06 Aug 12 peter 397
2809 06 Aug 12 peter 398        \note this function uses a biased variance estimate [1/n]
2809 06 Aug 12 peter 399        compared to kurtosis(utility::VectorBase&), which uses an
2809 06 Aug 12 peter 400        unbiased variance estimate [1/(n-1)] and may thus give
2809 06 Aug 12 peter 401        different results for small \c n.
2809 06 Aug 12 peter 402     */
2809 06 Aug 12 peter 403     double kurtosis(void) const
2809 06 Aug 12 peter 404     { return central_moment4()/std::pow(this->variance(),2) - 3.0; }
2809 06 Aug 12 peter 405
2809 06 Aug 12 peter 406     /**
2809 06 Aug 12 peter 407        \brief sum quadrupled centered
2809 06 Aug 12 peter 408
2809 06 Aug 12 peter 409        Calculated as \f$ \frac{1}{n} \sum (x_i-m)^4 \f$
2809 06 Aug 12 peter 410      */
2809 06 Aug 12 peter 411     double sum_x4_centered(void) const { return cm4_; }
2809 06 Aug 12 peter 412   protected:
2809 06 Aug 12 peter 413     /**
2809 06 Aug 12 peter 414        \brief sum quadrupled centered
2809 06 Aug 12 peter 415      */
2809 06 Aug 12 peter 416     double cm4_;
2809 06 Aug 12 peter 417
2809 06 Aug 12 peter 418     /**
2809 06 Aug 12 peter 419        add a set of \a n data points with mean \c mean, sum squared
2809 06 Aug 12 peter 420        centered \a cm2, sum cubed centered \a cm3, and sum quadrupled
2809 06 Aug 12 peter 421        centered \a cm4.
2809 06 Aug 12 peter 422     */
2809 06 Aug 12 peter 423     void add4(double mean, double cm2, double cm3, double cm4, long int n);
2809 06 Aug 12 peter 424
2809 06 Aug 12 peter 425     /**
2809 06 Aug 12 peter 426        add one data point with value \c delta + mean()
2809 06 Aug 12 peter 427      */
2809 06 Aug 12 peter 428     void add4(double delta);
2809 06 Aug 12 peter 429
2809 06 Aug 12 peter 430     /**
2809 06 Aug 12 peter 431        rescales \f$ cm4 = x^4 cm4 \f$ and calles
2809 06 Aug 12 peter 432        rescale2(double)
2809 06 Aug 12 peter 433      */
2809 06 Aug 12 peter 434     double rescale4(double x);
2809 06 Aug 12 peter 435   };
2809 06 Aug 12 peter 436
2809 06 Aug 12 peter 437   ///// template implementations
2809 06 Aug 12 peter 438
2809 06 Aug 12 peter 439   template<class Derived>
2809 06 Aug 12 peter 440   averager_base<Derived>::~averager_base(void) {}
2809 06 Aug 12 peter 441   template<class Derived>
2809 06 Aug 12 peter 442   averager_base2<Derived>::~averager_base2(void) {}
2809 06 Aug 12 peter 443   template<class Derived>
2809 06 Aug 12 peter 444   averager_base3<Derived>::~averager_base3(void) {}
2809 06 Aug 12 peter 445   template<class Derived>
2809 06 Aug 12 peter 446   averager_base4<Derived>::~averager_base4(void) {}
2809 06 Aug 12 peter 447
2809 06 Aug 12 peter 448   template<class Derived>
2809 06 Aug 12 peter 449   void averager_base<Derived>::add1(double x, long int n)
2809 06 Aug 12 peter 450   {
2809 06 Aug 12 peter 451     if (n == -this->n_) {
2809 06 Aug 12 peter 452       this->reset();
2809 06 Aug 12 peter 453       return;
2809 06 Aug 12 peter 454     }
2809 06 Aug 12 peter 455     this->mean_ = (this->n_*this->mean_ + n*x) / (this->n_ + n);
2809 06 Aug 12 peter 456     this->n_+=n;
2809 06 Aug 12 peter 457   }
2809 06 Aug 12 peter 458
2809 06 Aug 12 peter 459
2809 06 Aug 12 peter 460   template<class Derived>
2809 06 Aug 12 peter 461   void averager_base<Derived>::add1(double delta)
2809 06 Aug 12 peter 462   {
2809 06 Aug 12 peter 463     ++this->n_;
2809 06 Aug 12 peter 464     this->mean_ += delta * 1/(this->n_);
2809 06 Aug 12 peter 465   }
2809 06 Aug 12 peter 466
2809 06 Aug 12 peter 467
2809 06 Aug 12 peter 468   template<class Derived>
2809 06 Aug 12 peter 469   void averager_base2<Derived>::add2(double m, double cm2, long int n)
2809 06 Aug 12 peter 470   {
2809 06 Aug 12 peter 471     if (n == -this->n_) {
2809 06 Aug 12 peter 472       this->reset();
2809 06 Aug 12 peter 473       return;
2809 06 Aug 12 peter 474     }
2809 06 Aug 12 peter 475     const double delta = this->mean_ - m;
2809 06 Aug 12 peter 476     this->cm2_ += cm2 + this->n_*n*delta*delta/(this->n_+n);
2809 06 Aug 12 peter 477     this->add1(m, n);
2809 06 Aug 12 peter 478   }
2809 06 Aug 12 peter 479
2809 06 Aug 12 peter 480
2809 06 Aug 12 peter 481   template<class Derived>
2809 06 Aug 12 peter 482   void averager_base2<Derived>::add2(double delta)
2809 06 Aug 12 peter 483   {
2809 06 Aug 12 peter 484     this->cm2_ += this->n_*delta*delta/(this->n_+1);
2809 06 Aug 12 peter 485     this->add1(delta);
2809 06 Aug 12 peter 486   }
2809 06 Aug 12 peter 487
2809 06 Aug 12 peter 488
2809 06 Aug 12 peter 489   template<class Derived>
2809 06 Aug 12 peter 490   void averager_base3<Derived>::add3(double mean, double cm2, double cm3,
2809 06 Aug 12 peter 491                                      long int n)
2809 06 Aug 12 peter 492   {
2809 06 Aug 12 peter 493     if (n == -this->n_) {
2809 06 Aug 12 peter 494       this->reset();
2809 06 Aug 12 peter 495       return;
2809 06 Aug 12 peter 496     }
2809 06 Aug 12 peter 497     const double n_inverse = 1.0/(this->n_+n);
2809 06 Aug 12 peter 498     const double delta = mean - this->mean_;
2809 06 Aug 12 peter 499     const double r = n_inverse*delta;
2809 06 Aug 12 peter 500
2809 06 Aug 12 peter 501     this->cm3_ += cm3 + r*r*delta*this->n_*n*(this->n_-n) +
2809 06 Aug 12 peter 502       3*r*(this->n_*cm2 - n*this->cm2_);
2809 06 Aug 12 peter 503     this->add2(mean, cm2, n);
2809 06 Aug 12 peter 504   }
2809 06 Aug 12 peter 505
2809 06 Aug 12 peter 506
2809 06 Aug 12 peter 507   template<class Derived>
2809 06 Aug 12 peter 508   void averager_base3<Derived>::add3(double delta)
2809 06 Aug 12 peter 509   {
2809 06 Aug 12 peter 510     const double r = delta/(this->n_+1);
2809 06 Aug 12 peter 511
2809 06 Aug 12 peter 512     this->cm3_ += r*r*delta*this->n_*(this->n_-1) - 3*r*this->cm2_;
2809 06 Aug 12 peter 513     this->add2(delta);
2809 06 Aug 12 peter 514   }
2809 06 Aug 12 peter 515
2809 06 Aug 12 peter 516
2809 06 Aug 12 peter 517   template<class Derived>
2809 06 Aug 12 peter 518   void averager_base4<Derived>::add4(double mean, double cm2, double cm3,
2809 06 Aug 12 peter 519                                      double cm4, long int n)
2809 06 Aug 12 peter 520   {
2809 06 Aug 12 peter 521     if (n == -this->n_) {
2809 06 Aug 12 peter 522       this->reset();
2809 06 Aug 12 peter 523       return;
2809 06 Aug 12 peter 524     }
2809 06 Aug 12 peter 525     const double n_inverse = 1.0/(this->n_+n);
2809 06 Aug 12 peter 526     const double delta = mean - this->mean_;
2809 06 Aug 12 peter 527     const double r = n_inverse*delta;
2809 06 Aug 12 peter 528
2809 06 Aug 12 peter 529     cm4_ += cm4 +
2809 06 Aug 12 peter 530       std::pow(r,3)*delta*this->n_*n*(this->n_*this->n_-this->n_*n+n*n) +
2809 06 Aug 12 peter 531       6*r*r*(this->n_*this->n_*cm2+n*n*this->cm2_) +
2809 06 Aug 12 peter 532       4*r*(this->n_*cm3-n*this->cm3_);
2809 06 Aug 12 peter 533     this->add3(mean, cm2, cm3, n);
2809 06 Aug 12 peter 534   }
2809 06 Aug 12 peter 535
2809 06 Aug 12 peter 536
2809 06 Aug 12 peter 537   template<class Derived>
2809 06 Aug 12 peter 538   void averager_base4<Derived>::add4(double delta)
2809 06 Aug 12 peter 539   {
2809 06 Aug 12 peter 540     const double r = delta/(this->n_+1);
2809 06 Aug 12 peter 541
2809 06 Aug 12 peter 542     cm4_ +=
2809 06 Aug 12 peter 543       std::pow(r,3)*delta*this->n_*(this->n_*this->n_-this->n_+1) +
2809 06 Aug 12 peter 544       6*r*r*(this->cm2_) +
2809 06 Aug 12 peter 545       4*r*(0-this->cm3_);
2809 06 Aug 12 peter 546
2809 06 Aug 12 peter 547     this->add3(delta);
2809 06 Aug 12 peter 548   }
2809 06 Aug 12 peter 549
2809 06 Aug 12 peter 550
2809 06 Aug 12 peter 551   template<class Derived>
2809 06 Aug 12 peter 552   double averager_base<Derived>::rescale1(double x)
2809 06 Aug 12 peter 553   {
2809 06 Aug 12 peter 554     this->mean_ *= x;
2809 06 Aug 12 peter 555     return x;
2809 06 Aug 12 peter 556   }
2809 06 Aug 12 peter 557
2809 06 Aug 12 peter 558
2809 06 Aug 12 peter 559   template<class Derived>
2809 06 Aug 12 peter 560   double averager_base2<Derived>::rescale2(double x)
2809 06 Aug 12 peter 561   {
2809 06 Aug 12 peter 562     this->rescale1(x);
2809 06 Aug 12 peter 563     double factor = x*x;
2809 06 Aug 12 peter 564     this->cm2_ *= factor;
2809 06 Aug 12 peter 565     return factor;
2809 06 Aug 12 peter 566   }
2809 06 Aug 12 peter 567
2809 06 Aug 12 peter 568
2809 06 Aug 12 peter 569   template<class Derived>
2809 06 Aug 12 peter 570   double averager_base3<Derived>::rescale3(double x)
2809 06 Aug 12 peter 571   {
2809 06 Aug 12 peter 572     double factor = x*this->rescale2(x);
2809 06 Aug 12 peter 573     this->cm3_ *= factor;
2809 06 Aug 12 peter 574     return factor;
2809 06 Aug 12 peter 575   }
2809 06 Aug 12 peter 576
2809 06 Aug 12 peter 577
2809 06 Aug 12 peter 578   template<class Derived>
2809 06 Aug 12 peter 579   double averager_base4<Derived>::rescale4(double x)
2809 06 Aug 12 peter 580   {
2809 06 Aug 12 peter 581     double factor = x*this->rescale3(x);
2809 06 Aug 12 peter 582     this->cm3_ *= factor;
2809 06 Aug 12 peter 583     return factor;
2809 06 Aug 12 peter 584   }
2809 06 Aug 12 peter 585
2809 06 Aug 12 peter 586 }}}
2809 06 Aug 12 peter 587 #endif

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