yat/regression/LinearWeighted.cc

Code
Comments
Other
Rev Date Author Line
429 08 Dec 05 peter 1 // $Id$
429 08 Dec 05 peter 2
675 10 Oct 06 jari 3 /*
831 27 Mar 07 peter 4   Copyright (C) 2005 Peter Johansson
2119 12 Dec 09 peter 5   Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
4359 23 Aug 23 peter 6   Copyright (C) 2007 Peter Johansson
4359 23 Aug 23 peter 7   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 8   Copyright (C) 2012 Peter Johansson
429 08 Dec 05 peter 9
1437 25 Aug 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/yat
429 08 Dec 05 peter 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 "LinearWeighted.h"
682 11 Oct 06 jari 29 #include "yat/statistics/AveragerPairWeighted.h"
1120 21 Feb 08 peter 30 #include "yat/utility/Vector.h"
675 10 Oct 06 jari 31
730 06 Jan 07 peter 32 #include <cassert>
730 06 Jan 07 peter 33
429 08 Dec 05 peter 34 namespace theplu {
680 11 Oct 06 jari 35 namespace yat {
429 08 Dec 05 peter 36 namespace regression {
429 08 Dec 05 peter 37
703 18 Dec 06 jari 38   LinearWeighted::LinearWeighted(void)
703 18 Dec 06 jari 39     : OneDimensionalWeighted(), alpha_(0), alpha_var_(0), beta_(0),
730 06 Jan 07 peter 40       beta_var_(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   LinearWeighted::~LinearWeighted(void)
703 18 Dec 06 jari 45   {
703 18 Dec 06 jari 46   }
703 18 Dec 06 jari 47
718 26 Dec 06 jari 48   double LinearWeighted::alpha(void) const
718 26 Dec 06 jari 49   {
718 26 Dec 06 jari 50     return alpha_;
718 26 Dec 06 jari 51   }
718 26 Dec 06 jari 52
729 05 Jan 07 peter 53   double LinearWeighted::alpha_var(void) const
718 26 Dec 06 jari 54   {
730 06 Jan 07 peter 55     return alpha_var_;
718 26 Dec 06 jari 56   }
718 26 Dec 06 jari 57
718 26 Dec 06 jari 58   double LinearWeighted::beta(void) const
718 26 Dec 06 jari 59   {
718 26 Dec 06 jari 60     return beta_;
718 26 Dec 06 jari 61   }
718 26 Dec 06 jari 62
729 05 Jan 07 peter 63   double LinearWeighted::beta_var(void) const
718 26 Dec 06 jari 64   {
730 06 Jan 07 peter 65     return beta_var_;
718 26 Dec 06 jari 66   }
718 26 Dec 06 jari 67
1020 01 Feb 08 peter 68   void LinearWeighted::fit(const utility::VectorBase& x,
1020 01 Feb 08 peter 69                            const utility::VectorBase& y,
1020 01 Feb 08 peter 70                            const utility::VectorBase& w)
429 08 Dec 05 peter 71   {
730 06 Jan 07 peter 72     assert(x.size()==y.size());
730 06 Jan 07 peter 73     assert(x.size()==w.size());
730 06 Jan 07 peter 74
489 04 Jan 06 peter 75     // AveragerPairWeighted requires 2 weights but works only on the
489 04 Jan 06 peter 76     // product wx*wy, so we can send in w and a dummie to get what we
489 04 Jan 06 peter 77     // want.
702 26 Oct 06 peter 78     ap_.reset();
1120 21 Feb 08 peter 79     yat::utility::Vector dummy(x.size(), 1.0);
1043 06 Feb 08 peter 80     add(ap_, x.begin(), x.end(), y.begin(),dummy.begin(),w.begin());
489 04 Jan 06 peter 81
702 26 Oct 06 peter 82     alpha_ = m_y();
702 26 Oct 06 peter 83     beta_ = sxy()/sxx();
730 06 Jan 07 peter 84
730 06 Jan 07 peter 85     chisq_=0;
730 06 Jan 07 peter 86     for (size_t i=0; i<x.size(); ++i){
730 06 Jan 07 peter 87       double res = predict(x(i))-y(i);
730 06 Jan 07 peter 88       chisq_ += w(i)*res*res;
730 06 Jan 07 peter 89     }
730 06 Jan 07 peter 90
730 06 Jan 07 peter 91     alpha_var_ = s2()/ap_.y_averager().sum_w();
4200 19 Aug 22 peter 92     beta_var_ = s2()/sxx();
429 08 Dec 05 peter 93   }
429 08 Dec 05 peter 94
718 26 Dec 06 jari 95   double LinearWeighted::m_x(void) const
718 26 Dec 06 jari 96   {
718 26 Dec 06 jari 97     return ap_.x_averager().mean();
718 26 Dec 06 jari 98   }
718 26 Dec 06 jari 99
718 26 Dec 06 jari 100   double LinearWeighted::m_y(void) const
718 26 Dec 06 jari 101   {
718 26 Dec 06 jari 102     return ap_.y_averager().mean();
718 26 Dec 06 jari 103   }
718 26 Dec 06 jari 104
730 06 Jan 07 peter 105   double LinearWeighted::predict(const double x) const
4200 19 Aug 22 peter 106   {
4200 19 Aug 22 peter 107     return alpha_ + beta_ * (x-m_x());
718 26 Dec 06 jari 108   }
718 26 Dec 06 jari 109
4200 19 Aug 22 peter 110
718 26 Dec 06 jari 111   double LinearWeighted::s2(double w) const
718 26 Dec 06 jari 112   {
730 06 Jan 07 peter 113     return chisq_/(w*(ap_.y_averager().n()-2));
718 26 Dec 06 jari 114   }
718 26 Dec 06 jari 115
729 05 Jan 07 peter 116   double LinearWeighted::standard_error2(const double x) const
718 26 Dec 06 jari 117   {
730 06 Jan 07 peter 118     return alpha_var_ + beta_var_*(x-m_x())*(x-m_x());
718 26 Dec 06 jari 119   }
718 26 Dec 06 jari 120
730 06 Jan 07 peter 121
718 26 Dec 06 jari 122   double LinearWeighted::sxx(void) const
718 26 Dec 06 jari 123   {
718 26 Dec 06 jari 124     return ap_.x_averager().sum_xx_centered();
718 26 Dec 06 jari 125   }
718 26 Dec 06 jari 126
730 06 Jan 07 peter 127
718 26 Dec 06 jari 128   double LinearWeighted::sxy(void) const
718 26 Dec 06 jari 129   {
718 26 Dec 06 jari 130     return ap_.sum_xy_centered();
718 26 Dec 06 jari 131   }
718 26 Dec 06 jari 132
730 06 Jan 07 peter 133
718 26 Dec 06 jari 134   double LinearWeighted::syy(void) const
718 26 Dec 06 jari 135   {
718 26 Dec 06 jari 136     return ap_.y_averager().sum_xx_centered();
718 26 Dec 06 jari 137   }
718 26 Dec 06 jari 138
681 11 Oct 06 jari 139 }}} // of namespaces regression, yat, and theplu

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