yat/statistics/AveragerWeighted.cc

Code
Comments
Other
Rev Date Author Line
94 09 Jun 04 peter 1 // $Id$
94 09 Jun 04 peter 2
675 10 Oct 06 jari 3 /*
4359 23 Aug 23 peter 4   Copyright (C) 2004, 2005, 2006 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 5   Copyright (C) 2007 Peter Johansson
4359 23 Aug 23 peter 6   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 7   Copyright (C) 2011, 2012, 2014 Peter Johansson
94 09 Jun 04 peter 8
1437 25 Aug 08 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 10
675 10 Oct 06 jari 11   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 12   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 13   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 14   License, or (at your option) any later version.
675 10 Oct 06 jari 15
675 10 Oct 06 jari 16   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 19   General Public License for more details.
675 10 Oct 06 jari 20
675 10 Oct 06 jari 21   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 23 */
675 10 Oct 06 jari 24
2881 18 Nov 12 peter 25 #include <config.h>
2881 18 Nov 12 peter 26
680 11 Oct 06 jari 27 #include "AveragerWeighted.h"
675 10 Oct 06 jari 28
2673 03 Dec 11 peter 29 #include <limits>
2562 25 Sep 11 peter 30 #include <cmath>
2562 25 Sep 11 peter 31
94 09 Jun 04 peter 32 namespace theplu {
680 11 Oct 06 jari 33 namespace yat {
683 11 Oct 06 jari 34 namespace statistics {
94 09 Jun 04 peter 35
703 18 Dec 06 jari 36   AveragerWeighted::AveragerWeighted(void)
2562 25 Sep 11 peter 37     : mean_(0), m2_(0), w_(0), ww_(0), wwxx_(0), wxx_(0)
95 09 Jun 04 peter 38   {
703 18 Dec 06 jari 39   }
703 18 Dec 06 jari 40
2562 25 Sep 11 peter 41
703 18 Dec 06 jari 42   void AveragerWeighted::add(const double d, const double w)
703 18 Dec 06 jari 43   {
703 18 Dec 06 jari 44     if (!w)
703 18 Dec 06 jari 45       return;
2562 25 Sep 11 peter 46     double new_w = w_ + w;
2562 25 Sep 11 peter 47     double delta = d - mean_;
2562 25 Sep 11 peter 48     double R = delta * w / (new_w);
2562 25 Sep 11 peter 49     mean_ += R;
2562 25 Sep 11 peter 50     m2_ += w_ * delta * R;
2562 25 Sep 11 peter 51     w_ = new_w;
2562 25 Sep 11 peter 52     ww_ += w*w;
4200 19 Aug 22 peter 53     wwxx_ += w*w*d*d;
2562 25 Sep 11 peter 54     wxx_ += w*d*d;
703 18 Dec 06 jari 55   }
703 18 Dec 06 jari 56
2562 25 Sep 11 peter 57
718 26 Dec 06 jari 58   double AveragerWeighted::mean(void) const
718 26 Dec 06 jari 59   {
2562 25 Sep 11 peter 60     if (!w_) // no data, mean is undefined
2562 25 Sep 11 peter 61       return std::numeric_limits<double>::quiet_NaN();
2562 25 Sep 11 peter 62     return mean_;
718 26 Dec 06 jari 63   }
718 26 Dec 06 jari 64
2562 25 Sep 11 peter 65
718 26 Dec 06 jari 66   double AveragerWeighted::n(void) const
718 26 Dec 06 jari 67   {
3266 03 Jul 14 peter 68     if (sum_w()==0.0) // empty object, return 0
3266 03 Jul 14 peter 69       return 0.0;
718 26 Dec 06 jari 70     return sum_w()*sum_w()/sum_ww();
718 26 Dec 06 jari 71   }
718 26 Dec 06 jari 72
2562 25 Sep 11 peter 73
703 18 Dec 06 jari 74   void AveragerWeighted::rescale(double a)
703 18 Dec 06 jari 75   {
2562 25 Sep 11 peter 76     mean_ *= a;
2562 25 Sep 11 peter 77     double a2 = a*a;
2562 25 Sep 11 peter 78     m2_ *= a2;
2562 25 Sep 11 peter 79     wwxx_ *= a2;
2562 25 Sep 11 peter 80     wxx_*=a2;
703 18 Dec 06 jari 81   }
703 18 Dec 06 jari 82
2562 25 Sep 11 peter 83
703 18 Dec 06 jari 84   void AveragerWeighted::reset(void)
703 18 Dec 06 jari 85   {
2562 25 Sep 11 peter 86     AveragerWeighted tmp;
2562 25 Sep 11 peter 87     *this = tmp;
703 18 Dec 06 jari 88   }
703 18 Dec 06 jari 89
2562 25 Sep 11 peter 90
718 26 Dec 06 jari 91   double AveragerWeighted::std(void) const
718 26 Dec 06 jari 92   {
2562 25 Sep 11 peter 93     return std::sqrt(variance());
718 26 Dec 06 jari 94   }
718 26 Dec 06 jari 95
2562 25 Sep 11 peter 96
718 26 Dec 06 jari 97   double AveragerWeighted::standard_error(void) const
718 26 Dec 06 jari 98   {
2562 25 Sep 11 peter 99     return std::sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered());
718 26 Dec 06 jari 100   }
718 26 Dec 06 jari 101
2562 25 Sep 11 peter 102
718 26 Dec 06 jari 103   double AveragerWeighted::sum_w(void)  const
718 26 Dec 06 jari 104   {
2562 25 Sep 11 peter 105     return w_;
718 26 Dec 06 jari 106   }
718 26 Dec 06 jari 107
2562 25 Sep 11 peter 108
718 26 Dec 06 jari 109   double AveragerWeighted::sum_ww(void) const
718 26 Dec 06 jari 110   {
2562 25 Sep 11 peter 111     return ww_;
718 26 Dec 06 jari 112   }
718 26 Dec 06 jari 113
2562 25 Sep 11 peter 114
718 26 Dec 06 jari 115   double AveragerWeighted::sum_wwx(void) const
718 26 Dec 06 jari 116   {
2562 25 Sep 11 peter 117     return m2_ + mean_*mean_*w_;
718 26 Dec 06 jari 118   }
718 26 Dec 06 jari 119
2562 25 Sep 11 peter 120
718 26 Dec 06 jari 121   double AveragerWeighted::sum_wwxx(void) const
718 26 Dec 06 jari 122   {
2562 25 Sep 11 peter 123     return wwxx_;
718 26 Dec 06 jari 124   }
718 26 Dec 06 jari 125
2562 25 Sep 11 peter 126
718 26 Dec 06 jari 127   double AveragerWeighted::sum_wx(void) const
718 26 Dec 06 jari 128   {
2562 25 Sep 11 peter 129     if (!w_)
2562 25 Sep 11 peter 130       return 0.0;
2562 25 Sep 11 peter 131     return mean_ * w_;
718 26 Dec 06 jari 132   }
718 26 Dec 06 jari 133
2562 25 Sep 11 peter 134
718 26 Dec 06 jari 135   double AveragerWeighted::sum_wxx(void) const
718 26 Dec 06 jari 136   {
718 26 Dec 06 jari 137     return wxx_;
718 26 Dec 06 jari 138   }
718 26 Dec 06 jari 139
2562 25 Sep 11 peter 140
718 26 Dec 06 jari 141   double AveragerWeighted::sum_xx_centered(void) const
718 26 Dec 06 jari 142   {
2562 25 Sep 11 peter 143     return m2_;
718 26 Dec 06 jari 144   }
718 26 Dec 06 jari 145
2562 25 Sep 11 peter 146
718 26 Dec 06 jari 147   double AveragerWeighted::variance(const double m) const
718 26 Dec 06 jari 148   {
718 26 Dec 06 jari 149     return (sum_wxx()-2*m*sum_wx())/sum_w()+m*m;
718 26 Dec 06 jari 150   }
718 26 Dec 06 jari 151
2562 25 Sep 11 peter 152
718 26 Dec 06 jari 153   double AveragerWeighted::variance(void) const
718 26 Dec 06 jari 154   {
718 26 Dec 06 jari 155     return sum_xx_centered()/sum_w();
718 26 Dec 06 jari 156   }
718 26 Dec 06 jari 157
718 26 Dec 06 jari 158
4200 19 Aug 22 peter 159   const AveragerWeighted&
2562 25 Sep 11 peter 160   AveragerWeighted::operator+=(const AveragerWeighted& a)
718 26 Dec 06 jari 161   {
2563 25 Sep 11 peter 162     if (!a.w_)
2563 25 Sep 11 peter 163       return *this;
2562 25 Sep 11 peter 164     double delta = mean() - a.mean();
2798 27 Jul 12 peter 165     mean_ = (sum_wx() + a.sum_wx()) / (sum_w() + a.sum_w());
2562 25 Sep 11 peter 166     m2_ += a.m2_ + sum_w()*a.sum_w()*delta*delta/(sum_w()+a.sum_w());
2562 25 Sep 11 peter 167     w_ += a.w_;
2562 25 Sep 11 peter 168     ww_ += a.ww_;
2562 25 Sep 11 peter 169     wwxx_ += a.wwxx_;
2562 25 Sep 11 peter 170     wxx_ += a.wxx_;
397 13 Oct 05 jari 171     return *this;
95 09 Jun 04 peter 172   }
95 09 Jun 04 peter 173
683 11 Oct 06 jari 174 }}} // of namespace statistics, yat, and theplu

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