source: trunk/yat/statistics/AveragerWeighted.cc

Last change on this file was 4245, checked in by Peter, 2 weeks ago

GCC -Wdeprecated-copy complains when some but not all of copy ctor, dtor and operator= are user-implemented. When the implementation is identical with the one generated by the compiler, there is no reason to keep the code. DataWeightProxy? is a different story (see #986)

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