source: trunk/yat/statistics/AveragerWeighted.cc

Last change on this file was 3330, checked in by Peter, 7 years ago

update copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.3 KB
Line 
1// $Id: AveragerWeighted.cc 3330 2014-10-14 08:03:25Z peter $
2
3/*
4  Copyright (C) 2004, 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2011, 2012, 2014 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  AveragerWeighted::AveragerWeighted(const AveragerWeighted& a)
41    : mean_(a.mean_), m2_(a.m2_), w_(a.w_),
42      ww_(a.ww_), wwxx_(a.wwxx_), wxx_(a.wxx_)
43  {
44  }
45
46
47  void AveragerWeighted::add(const double d, const double w)
48  {
49    if (!w)
50      return;
51    double new_w = w_ + w;
52    double delta = d - mean_;
53    double R = delta * w / (new_w);
54    mean_ += R;
55    m2_ += w_ * delta * R;
56    w_ = new_w;
57    ww_ += w*w;
58    wwxx_ += w*w*d*d; 
59    wxx_ += w*d*d;
60  }
61
62
63  double AveragerWeighted::mean(void) const
64  {
65    if (!w_) // no data, mean is undefined
66      return std::numeric_limits<double>::quiet_NaN();
67    return mean_;
68  }
69
70
71  double AveragerWeighted::n(void) const
72  {
73    if (sum_w()==0.0) // empty object, return 0
74      return 0.0;
75    return sum_w()*sum_w()/sum_ww();
76  }
77
78
79  void AveragerWeighted::rescale(double a)
80  {
81    mean_ *= a;
82    double a2 = a*a;
83    m2_ *= a2;
84    wwxx_ *= a2;
85    wxx_*=a2;
86  }
87
88
89  void AveragerWeighted::reset(void)
90  {
91    AveragerWeighted tmp;
92    *this = tmp;
93  }
94
95
96  double AveragerWeighted::std(void) const
97  {
98    return std::sqrt(variance());
99  }
100
101
102  double AveragerWeighted::standard_error(void) const
103  {
104    return std::sqrt(sum_ww()/(sum_w()*sum_w()*sum_w()) * sum_xx_centered());
105  }
106
107
108  double AveragerWeighted::sum_w(void)  const
109  {
110    return w_;
111  }
112
113
114  double AveragerWeighted::sum_ww(void) const
115  {
116    return ww_;
117  }
118
119
120  double AveragerWeighted::sum_wwx(void) const
121  {
122    return m2_ + mean_*mean_*w_;
123  }
124
125
126  double AveragerWeighted::sum_wwxx(void) const
127  {
128    return wwxx_;
129  }
130
131
132  double AveragerWeighted::sum_wx(void) const
133  {
134    if (!w_)
135      return 0.0;
136    return mean_ * w_;
137  }
138
139
140  double AveragerWeighted::sum_wxx(void) const
141  {
142    return wxx_;
143  }
144
145
146  double AveragerWeighted::sum_xx_centered(void) const
147  {
148    return m2_;
149  }
150
151
152  double AveragerWeighted::variance(const double m) const
153  {
154    return (sum_wxx()-2*m*sum_wx())/sum_w()+m*m;
155  }
156
157
158  double AveragerWeighted::variance(void) const
159  {
160    return sum_xx_centered()/sum_w();
161  }
162
163
164  const AveragerWeighted& 
165  AveragerWeighted::operator+=(const AveragerWeighted& a)
166  {
167    if (!a.w_)
168      return *this;
169    double delta = mean() - a.mean();
170    mean_ = (sum_wx() + a.sum_wx()) / (sum_w() + a.sum_w());
171    m2_ += a.m2_ + sum_w()*a.sum_w()*delta*delta/(sum_w()+a.sum_w());
172    w_ += a.w_;
173    ww_ += a.ww_;
174    wwxx_ += a.wwxx_;
175    wxx_ += a.wxx_;
176    return *this;
177  }
178
179}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.