yat  0.14.5pre
averager_base.h
1 #ifndef _theplu_yat_statistics_average_base
2 #define _theplu_yat_statistics_average_base
3 
4 // $Id: averager_base.h 2809 2012-08-06 13:44:47Z peter $
5 
6 /*
7  Copyright (C) 2012 Peter Johansson
8 
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10 
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15 
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with yat. If not, see <http://www.gnu.org/licenses/>.
23 */
24 
25 #include <boost/concept_check.hpp>
26 
27 #include <cmath>
28 #include <limits>
29 
30 namespace theplu {
31 namespace yat {
32 namespace statistics {
33 
68  template<class Derived>
70  {
71  public:
75  averager_base(void) : n_(0), mean_(0)
76  {
77  BOOST_CONCEPT_ASSERT((boost::DefaultConstructible<Derived>));
78  }
79 
83  virtual ~averager_base(void)=0;
84 
90  void add(double x, long n=1) {static_cast<Derived*>(this)->add_impl(x, n);}
91 
97  double mean(void) const
98  { return n_ ? mean_ : std::numeric_limits<double>::quiet_NaN(); }
99 
105  long n(void) const { return n_; }
106 
112  void rescale(double a) { static_cast<Derived*>(this)->rescale_impl(a); }
113 
119  void reset(void)
120  {
121  Derived tmp;
122  *static_cast<Derived*>(this) = tmp;
123  }
124 
128  double sum_x(void) const { return n_*mean_; }
129 
130  protected:
134  long int n_;
135 
139  double mean_;
140 
144  void add1(double x, long int n);
148  void add1(double delta);
149 
153  double rescale1(double x);
154  };
155 
156 
169  template<class Derived>
170  class averager_base2 : public averager_base<Derived>
171  {
172  public:
176  averager_base2(void) : cm2_(0) {}
177 
181  virtual ~averager_base2(void)=0;
182 
191  double cv(void) const { return std() / this->mean(); }
192 
197  double standard_error(void) const
198  { return std::sqrt(variance()/this->n()); }
199 
206  double std(void) const { return std::sqrt(variance()); }
207 
214  double std(double m) const { return std::sqrt(variance(m)); }
215 
219  double sum_xx(void) const
220  { return cm2_ + this->n_*this->mean_*this->mean_; }
221 
225  double sum_xx_centered(void) const { return cm2_; }
226 
235  double variance(double m) const
236  { return sum_xx()/this->n() + m*m - 2*m*this->mean();}
237 
246  double variance(void) const { return sum_xx_centered()/this->n_; }
247 
256  double variance_unbiased(void) const
257  {
258  return (this->n_>1) ? sum_xx_centered()/(this->n_-1) :
259  std::numeric_limits<double>::quiet_NaN();
260  }
261 
262  protected:
266  double cm2_;
267 
272  void add2(double mean, double cm2, long int n);
273 
277  void add2(double delta);
278 
282  double rescale2(double x);
283  private:
284  };
285 
286 
295  template<class Derived>
296  class averager_base3 : public averager_base2<Derived>
297  {
298  public:
302  averager_base3(void) : averager_base2<Derived>(), cm3_(0) {}
303 
307  virtual ~averager_base3(void)=0;
308 
315  double central_moment3(void) const { return cm3_/this->n_; }
316 
328  double skewness(void) const
329  { return central_moment3()/std::pow(this->std(),3); }
330 
336  double sum_x3_centered(void) const { return cm3_; }
337 
338  protected:
342  double cm3_;
343 
348  void add3(double mean, double cm2, double cm3, long int n);
349 
353  void add3(double delta);
354 
358  double rescale3(double x);
359  };
360 
369  template<class Derived>
370  class averager_base4 : public averager_base3<Derived>
371  {
372  public:
376  averager_base4(void) : averager_base3<Derived>(), cm4_(0) {}
377 
381  virtual ~averager_base4(void)=0;
382 
389  double central_moment4(void) const { return cm4_/this->n_; }
390 
403  double kurtosis(void) const
404  { return central_moment4()/std::pow(this->variance(),2) - 3.0; }
405 
411  double sum_x4_centered(void) const { return cm4_; }
412  protected:
416  double cm4_;
417 
423  void add4(double mean, double cm2, double cm3, double cm4, long int n);
424 
428  void add4(double delta);
429 
434  double rescale4(double x);
435  };
436 
438 
439  template<class Derived>
441  template<class Derived>
443  template<class Derived>
445  template<class Derived>
447 
448  template<class Derived>
449  void averager_base<Derived>::add1(double x, long int n)
450  {
451  if (n == -this->n_) {
452  this->reset();
453  return;
454  }
455  this->mean_ = (this->n_*this->mean_ + n*x) / (this->n_ + n);
456  this->n_+=n;
457  }
458 
459 
460  template<class Derived>
461  void averager_base<Derived>::add1(double delta)
462  {
463  ++this->n_;
464  this->mean_ += delta * 1/(this->n_);
465  }
466 
467 
468  template<class Derived>
469  void averager_base2<Derived>::add2(double m, double cm2, long int n)
470  {
471  if (n == -this->n_) {
472  this->reset();
473  return;
474  }
475  const double delta = this->mean_ - m;
476  this->cm2_ += cm2 + this->n_*n*delta*delta/(this->n_+n);
477  this->add1(m, n);
478  }
479 
480 
481  template<class Derived>
482  void averager_base2<Derived>::add2(double delta)
483  {
484  this->cm2_ += this->n_*delta*delta/(this->n_+1);
485  this->add1(delta);
486  }
487 
488 
489  template<class Derived>
490  void averager_base3<Derived>::add3(double mean, double cm2, double cm3,
491  long int n)
492  {
493  if (n == -this->n_) {
494  this->reset();
495  return;
496  }
497  const double n_inverse = 1.0/(this->n_+n);
498  const double delta = mean - this->mean_;
499  const double r = n_inverse*delta;
500 
501  this->cm3_ += cm3 + r*r*delta*this->n_*n*(this->n_-n) +
502  3*r*(this->n_*cm2 - n*this->cm2_);
503  this->add2(mean, cm2, n);
504  }
505 
506 
507  template<class Derived>
508  void averager_base3<Derived>::add3(double delta)
509  {
510  const double r = delta/(this->n_+1);
511 
512  this->cm3_ += r*r*delta*this->n_*(this->n_-1) - 3*r*this->cm2_;
513  this->add2(delta);
514  }
515 
516 
517  template<class Derived>
518  void averager_base4<Derived>::add4(double mean, double cm2, double cm3,
519  double cm4, long int n)
520  {
521  if (n == -this->n_) {
522  this->reset();
523  return;
524  }
525  const double n_inverse = 1.0/(this->n_+n);
526  const double delta = mean - this->mean_;
527  const double r = n_inverse*delta;
528 
529  cm4_ += cm4 +
530  std::pow(r,3)*delta*this->n_*n*(this->n_*this->n_-this->n_*n+n*n) +
531  6*r*r*(this->n_*this->n_*cm2+n*n*this->cm2_) +
532  4*r*(this->n_*cm3-n*this->cm3_);
533  this->add3(mean, cm2, cm3, n);
534  }
535 
536 
537  template<class Derived>
538  void averager_base4<Derived>::add4(double delta)
539  {
540  const double r = delta/(this->n_+1);
541 
542  cm4_ +=
543  std::pow(r,3)*delta*this->n_*(this->n_*this->n_-this->n_+1) +
544  6*r*r*(this->cm2_) +
545  4*r*(0-this->cm3_);
546 
547  this->add3(delta);
548  }
549 
550 
551  template<class Derived>
553  {
554  this->mean_ *= x;
555  return x;
556  }
557 
558 
559  template<class Derived>
561  {
562  this->rescale1(x);
563  double factor = x*x;
564  this->cm2_ *= factor;
565  return factor;
566  }
567 
568 
569  template<class Derived>
571  {
572  double factor = x*this->rescale2(x);
573  this->cm3_ *= factor;
574  return factor;
575  }
576 
577 
578  template<class Derived>
580  {
581  double factor = x*this->rescale3(x);
582  this->cm3_ *= factor;
583  return factor;
584  }
585 
586 }}}
587 #endif
double rescale1(double x)
Definition: averager_base.h:552
double cm3_
Definition: averager_base.h:342
Definition: averager_base.h:370
averager_base2(void)
Constructor.
Definition: averager_base.h:176
double std(void) const
The standard deviation is defined as the square root of the variance.
Definition: averager_base.h:206
averager_base4(void)
Constructor.
Definition: averager_base.h:376
double cv(void) const
Coeffient of variation.
Definition: averager_base.h:191
virtual ~averager_base4(void)=0
Destructor.
Definition: averager_base.h:446
long int n_
Definition: averager_base.h:134
void add4(double mean, double cm2, double cm3, double cm4, long int n)
Definition: averager_base.h:518
double mean_
mean
Definition: averager_base.h:139
double rescale3(double x)
Definition: averager_base.h:570
Definition: averager_base.h:296
double central_moment3(void) const
Third central moment.
Definition: averager_base.h:315
void rescale(double a)
Rescales the object.
Definition: averager_base.h:112
Base class for averager classes.
Definition: averager_base.h:69
virtual ~averager_base3(void)=0
Destructor.
Definition: averager_base.h:444
void add2(double mean, double cm2, long int n)
Definition: averager_base.h:469
void reset(void)
Reset object.
Definition: averager_base.h:119
double variance(void) const
The estimated variance.
Definition: averager_base.h:246
double sum_x3_centered(void) const
Definition: averager_base.h:336
long n(void) const
number of data points
Definition: averager_base.h:105
Base class for averagers calculating mean and variance.
Definition: averager_base.h:170
void add(double x, long n=1)
add a data point
Definition: averager_base.h:90
double central_moment4(void) const
Fourth central moment.
Definition: averager_base.h:389
double rescale2(double x)
Definition: averager_base.h:560
double variance(double m) const
The variance with known mean.
Definition: averager_base.h:235
double skewness(void) const
Skewness.
Definition: averager_base.h:328
virtual ~averager_base(void)=0
Destructor.
Definition: averager_base.h:440
double sum_x4_centered(void) const
sum quadrupled centered
Definition: averager_base.h:411
void add1(double x, long int n)
Definition: averager_base.h:449
double sum_xx(void) const
Definition: averager_base.h:219
double mean(void) const
mean
Definition: averager_base.h:97
double rescale4(double x)
Definition: averager_base.h:579
double std(double m) const
The standard deviation is defined as the square root of the variance.
Definition: averager_base.h:214
double sum_xx_centered(void) const
Definition: averager_base.h:225
void add3(double mean, double cm2, double cm3, long int n)
Definition: averager_base.h:490
double cm2_
Definition: averager_base.h:266
double variance_unbiased(void) const
Definition: averager_base.h:256
averager_base3(void)
Constructor.
Definition: averager_base.h:302
virtual ~averager_base2(void)=0
Destructor.
Definition: averager_base.h:442
double cm4_
sum quadrupled centered
Definition: averager_base.h:416
double standard_error(void) const
Definition: averager_base.h:197
double kurtosis(void) const
Kurtosis.
Definition: averager_base.h:403
double sum_x(void) const
Definition: averager_base.h:128
averager_base(void)
Constructor.
Definition: averager_base.h:75

Generated on Tue Sep 26 2017 02:33:29 for yat by  doxygen 1.8.5