yat/statistics/AveragerPairWeighted.cc

Code
Comments
Other
Rev Date Author Line
474 22 Dec 05 peter 1 // $Id$
474 22 Dec 05 peter 2
675 10 Oct 06 jari 3 /*
2119 12 Dec 09 peter 4   Copyright (C) 2005 Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 5   Copyright (C) 2006 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 6   Copyright (C) 2007 Peter Johansson, Markus Ringnér
2119 12 Dec 09 peter 7   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 8   Copyright (C) 2011, 2012 Peter Johansson
474 22 Dec 05 peter 9
1437 25 Aug 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 11
675 10 Oct 06 jari 12   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 13   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 14   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 15   License, or (at your option) any later version.
675 10 Oct 06 jari 16
675 10 Oct 06 jari 17   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 18   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 20   General Public License for more details.
675 10 Oct 06 jari 21
675 10 Oct 06 jari 22   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 23   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 24 */
675 10 Oct 06 jari 25
2881 18 Nov 12 peter 26 #include <config.h>
2881 18 Nov 12 peter 27
680 11 Oct 06 jari 28 #include "AveragerPairWeighted.h"
680 11 Oct 06 jari 29 #include "AveragerPair.h"
680 11 Oct 06 jari 30 #include "Averager.h"
675 10 Oct 06 jari 31
781 05 Mar 07 peter 32 #include <cassert>
781 05 Mar 07 peter 33
474 22 Dec 05 peter 34 namespace theplu {
680 11 Oct 06 jari 35 namespace yat {
680 11 Oct 06 jari 36 namespace statistics {
474 22 Dec 05 peter 37
703 18 Dec 06 jari 38
703 18 Dec 06 jari 39   AveragerPairWeighted::AveragerPairWeighted(void)
2568 26 Sep 11 peter 40     : wxy_centered_(0)
703 18 Dec 06 jari 41   {
703 18 Dec 06 jari 42   }
703 18 Dec 06 jari 43
703 18 Dec 06 jari 44
4200 19 Aug 22 peter 45   void  AveragerPairWeighted::add(const double x, const double y,
474 22 Dec 05 peter 46                                   const double wx, const double wy)
4200 19 Aug 22 peter 47   {
476 22 Dec 05 markus 48     if(wx==0.0 || wy==0.0) {
476 22 Dec 05 markus 49       return;
476 22 Dec 05 markus 50     }
916 30 Sep 07 peter 51     assert(!std::isnan(x) && "x is nan");
916 30 Sep 07 peter 52     assert(!std::isnan(y) && "y is nan");
916 30 Sep 07 peter 53     assert(!std::isnan(wx) && "wx is nan");
916 30 Sep 07 peter 54     assert(!std::isnan(wy) && "wy is nan");
700 26 Oct 06 peter 55     double w=wx*wy;
2568 26 Sep 11 peter 56     xy_add(x, y, 0, w);
474 22 Dec 05 peter 57     x_.add(x,w);
2567 26 Sep 11 peter 58     y_.add(y,w);
474 22 Dec 05 peter 59   }
474 22 Dec 05 peter 60
627 05 Sep 06 peter 61
718 26 Dec 06 jari 62   double AveragerPairWeighted::correlation(void) const
718 26 Dec 06 jari 63   {
2569 26 Sep 11 peter 64     return sum_xy_centered() / sqrt(x_.sum_xx_centered()*y_.sum_xx_centered());
718 26 Dec 06 jari 65   }
718 26 Dec 06 jari 66
718 26 Dec 06 jari 67
718 26 Dec 06 jari 68   double AveragerPairWeighted::covariance(void) const
718 26 Dec 06 jari 69   {
718 26 Dec 06 jari 70     return sum_xy_centered()/sum_w();
718 26 Dec 06 jari 71   }
718 26 Dec 06 jari 72
718 26 Dec 06 jari 73
772 25 Feb 07 peter 74   double AveragerPairWeighted::msd(void) const
772 25 Feb 07 peter 75   {
2568 26 Sep 11 peter 76     return (x_.sum_wxx()+y_.sum_wxx()-2*sum_xy())/sum_w();
772 25 Feb 07 peter 77   }
772 25 Feb 07 peter 78
772 25 Feb 07 peter 79
899 26 Sep 07 peter 80   double AveragerPairWeighted::n(void) const
703 18 Dec 06 jari 81   {
899 26 Sep 07 peter 82     return x_.n();
703 18 Dec 06 jari 83   }
703 18 Dec 06 jari 84
899 26 Sep 07 peter 85
899 26 Sep 07 peter 86   void AveragerPairWeighted::reset(void)
890 25 Sep 07 markus 87   {
2568 26 Sep 11 peter 88     x_.reset(); y_.reset(); wxy_centered_=0;
890 25 Sep 07 markus 89   }
890 25 Sep 07 markus 90
899 26 Sep 07 peter 91
718 26 Dec 06 jari 92   double AveragerPairWeighted::sum_w(void) const
718 26 Dec 06 jari 93   {
2567 26 Sep 11 peter 94     return x_.sum_w();
718 26 Dec 06 jari 95   }
718 26 Dec 06 jari 96
718 26 Dec 06 jari 97
718 26 Dec 06 jari 98   double AveragerPairWeighted::sum_xy(void) const
718 26 Dec 06 jari 99   {
2568 26 Sep 11 peter 100     return sum_xy_centered() + x_.sum_wx()*y_.mean();
718 26 Dec 06 jari 101   }
718 26 Dec 06 jari 102
718 26 Dec 06 jari 103
718 26 Dec 06 jari 104   double AveragerPairWeighted::sum_xy_centered(void) const
718 26 Dec 06 jari 105   {
2568 26 Sep 11 peter 106     return wxy_centered_;
718 26 Dec 06 jari 107   }
718 26 Dec 06 jari 108
718 26 Dec 06 jari 109
4200 19 Aug 22 peter 110   void AveragerPairWeighted::xy_add(double mx, double my, double wxy_centered,
2568 26 Sep 11 peter 111                                     double w)
2568 26 Sep 11 peter 112   {
2568 26 Sep 11 peter 113     if (!w)
2568 26 Sep 11 peter 114       return;
2568 26 Sep 11 peter 115     /*
4200 19 Aug 22 peter 116       For derivation of update formula refer to
2568 26 Sep 11 peter 117       http://www.janinebennett.org/index_files/ParallelStatisticsAlgorithms.pdf
2568 26 Sep 11 peter 118     */
2568 26 Sep 11 peter 119     wxy_centered_ += wxy_centered;
2568 26 Sep 11 peter 120     if (sum_w())
4200 19 Aug 22 peter 121       wxy_centered_ +=
2568 26 Sep 11 peter 122         w*sum_w() / (w+sum_w()) * (mx-x_.mean()) * (my-y_.mean());
2568 26 Sep 11 peter 123   }
2568 26 Sep 11 peter 124
2568 26 Sep 11 peter 125
718 26 Dec 06 jari 126   const AveragerWeighted& AveragerPairWeighted::x_averager(void) const
718 26 Dec 06 jari 127   {
718 26 Dec 06 jari 128     return x_;
718 26 Dec 06 jari 129   }
718 26 Dec 06 jari 130
718 26 Dec 06 jari 131
718 26 Dec 06 jari 132   const AveragerWeighted& AveragerPairWeighted::y_averager(void) const
718 26 Dec 06 jari 133   {
718 26 Dec 06 jari 134     return y_;
718 26 Dec 06 jari 135   }
718 26 Dec 06 jari 136
4200 19 Aug 22 peter 137   const AveragerPairWeighted&
2570 26 Sep 11 peter 138   AveragerPairWeighted::operator+=(const AveragerPairWeighted& rhs)
2570 26 Sep 11 peter 139   {
4200 19 Aug 22 peter 140     xy_add(rhs.x_averager().mean(), rhs.y_averager().mean(),
2570 26 Sep 11 peter 141            rhs.sum_xy_centered(), rhs.sum_w());
2570 26 Sep 11 peter 142     x_ += rhs.x_averager();
2570 26 Sep 11 peter 143     y_ += rhs.y_averager();
2570 26 Sep 11 peter 144     return *this;
2570 26 Sep 11 peter 145   }
2570 26 Sep 11 peter 146
2570 26 Sep 11 peter 147
680 11 Oct 06 jari 148 }}} // of namespace statistics, yat, and theplu

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