yat/regression/PolynomialWeighted.cc

Code
Comments
Other
Rev Date Author Line
586 19 Jun 06 peter 1 // $Id$
586 19 Jun 06 peter 2
675 10 Oct 06 jari 3 /*
4359 23 Aug 23 peter 4   Copyright (C) 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) 2012 Peter Johansson
586 19 Jun 06 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 "PolynomialWeighted.h"
1121 22 Feb 08 peter 28 #include "yat/utility/Matrix.h"
1120 21 Feb 08 peter 29 #include "yat/utility/Vector.h"
675 10 Oct 06 jari 30
586 19 Jun 06 peter 31 #include <cassert>
586 19 Jun 06 peter 32
586 19 Jun 06 peter 33 namespace theplu {
680 11 Oct 06 jari 34 namespace yat {
586 19 Jun 06 peter 35 namespace regression {
586 19 Jun 06 peter 36
703 18 Dec 06 jari 37   PolynomialWeighted::PolynomialWeighted(size_t power)
740 12 Jan 07 peter 38     : OneDimensionalWeighted(), power_(power)
703 18 Dec 06 jari 39   {
703 18 Dec 06 jari 40   }
703 18 Dec 06 jari 41
703 18 Dec 06 jari 42   PolynomialWeighted::~PolynomialWeighted(void)
703 18 Dec 06 jari 43   {
703 18 Dec 06 jari 44   }
703 18 Dec 06 jari 45
1020 01 Feb 08 peter 46   void PolynomialWeighted::fit(const utility::VectorBase& x,
1020 01 Feb 08 peter 47                                const utility::VectorBase& y,
1020 01 Feb 08 peter 48                                const utility::VectorBase& w)
586 19 Jun 06 peter 49   {
586 19 Jun 06 peter 50     assert(x.size()==y.size());
586 19 Jun 06 peter 51     assert(y.size()==w.size());
741 13 Jan 07 peter 52     ap_.reset();
741 13 Jan 07 peter 53     // AveragerPairWeighted requires 2 weights but works only on the
741 13 Jan 07 peter 54     // product wx*wy, so we can send in w and a dummie to get what we
741 13 Jan 07 peter 55     // want.
1120 21 Feb 08 peter 56     utility::Vector dummy(x.size(), 1.0);
1043 06 Feb 08 peter 57     add(ap_,x.begin(), x.end(),y.begin(),dummy.begin(),w.begin());
4293 04 Feb 23 peter 58     utility::Matrix X(x.size(), power_+1, 1);
586 19 Jun 06 peter 59     for (size_t i=0; i<X.rows(); ++i)
1271 09 Apr 08 peter 60       for (size_t j=1; j<X.columns(); ++j)
586 19 Jun 06 peter 61         X(i,j)=X(i,j-1)*x(i);
586 19 Jun 06 peter 62     md_.fit(X,y,w);
741 13 Jan 07 peter 63     chisq_=md_.chisq();
586 19 Jun 06 peter 64   }
586 19 Jun 06 peter 65
740 12 Jan 07 peter 66
1120 21 Feb 08 peter 67   const utility::Vector& PolynomialWeighted::fit_parameters(void) const
718 26 Dec 06 jari 68   {
4200 19 Aug 22 peter 69     return md_.fit_parameters();
718 26 Dec 06 jari 70   }
718 26 Dec 06 jari 71
740 12 Jan 07 peter 72
742 13 Jan 07 peter 73   double PolynomialWeighted::s2(const double w) const
740 12 Jan 07 peter 74   {
740 12 Jan 07 peter 75     return md_.s2(w);
740 12 Jan 07 peter 76   }
740 12 Jan 07 peter 77
4293 04 Feb 23 peter 78
586 19 Jun 06 peter 79   double PolynomialWeighted::predict(const double x) const
586 19 Jun 06 peter 80   {
1120 21 Feb 08 peter 81     utility::Vector vec(power_+1,1);
586 19 Jun 06 peter 82     for (size_t i=1; i<=power_; ++i)
586 19 Jun 06 peter 83       vec(i) = vec(i-1)*x;
586 19 Jun 06 peter 84     return md_.predict(vec);
586 19 Jun 06 peter 85   }
586 19 Jun 06 peter 86
4293 04 Feb 23 peter 87
729 05 Jan 07 peter 88   double PolynomialWeighted::standard_error2(const double x) const
586 19 Jun 06 peter 89   {
1120 21 Feb 08 peter 90     utility::Vector vec(power_+1,1);
586 19 Jun 06 peter 91     for (size_t i=1; i<=power_; ++i)
586 19 Jun 06 peter 92       vec(i) = vec(i-1)*x;
731 06 Jan 07 peter 93     return md_.standard_error2(vec);
586 19 Jun 06 peter 94   }
586 19 Jun 06 peter 95
681 11 Oct 06 jari 96 }}} // of namespaces regression, yat, and theplu