yat/regression/Polynomial.cc

Code
Comments
Other
Rev Date Author Line
383 12 Aug 05 jari 1 // $Id$
383 12 Aug 05 jari 2
675 10 Oct 06 jari 3 /*
2119 12 Dec 09 peter 4   Copyright (C) 2005 Jari Häkkinen
4359 23 Aug 23 peter 5   Copyright (C) 2006 Jari Häkkinen, Peter Johansson
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
383 12 Aug 05 jari 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 "Polynomial.h"
1121 22 Feb 08 peter 29 #include "yat/utility/Matrix.h"
1019 01 Feb 08 peter 30 #include "yat/utility/VectorBase.h"
675 10 Oct 06 jari 31
383 12 Aug 05 jari 32 namespace theplu {
680 11 Oct 06 jari 33 namespace yat {
383 12 Aug 05 jari 34 namespace regression {
383 12 Aug 05 jari 35
703 18 Dec 06 jari 36   Polynomial::Polynomial(size_t power)
713 21 Dec 06 peter 37     : OneDimensional(), power_(power)
703 18 Dec 06 jari 38   {
703 18 Dec 06 jari 39   }
703 18 Dec 06 jari 40
713 21 Dec 06 peter 41
703 18 Dec 06 jari 42   Polynomial::~Polynomial(void)
703 18 Dec 06 jari 43   {
703 18 Dec 06 jari 44   }
703 18 Dec 06 jari 45
713 21 Dec 06 peter 46
1121 22 Feb 08 peter 47   const utility::Matrix& Polynomial::covariance(void) const
726 04 Jan 07 peter 48   {
726 04 Jan 07 peter 49     return md_.covariance();
726 04 Jan 07 peter 50   }
726 04 Jan 07 peter 51
726 04 Jan 07 peter 52
4200 19 Aug 22 peter 53   void Polynomial::fit(const utility::VectorBase& x,
1019 01 Feb 08 peter 54                        const utility::VectorBase& y)
383 12 Aug 05 jari 55   {
1043 06 Feb 08 peter 56     add(ap_, x.begin(), x.end(), y.begin());
2533 31 Jul 11 peter 57     utility::Matrix X(x.size(),power_+1,1);
386 14 Aug 05 jari 58     for (size_t i=0; i<X.rows(); ++i)
1271 09 Apr 08 peter 59       for (size_t j=1; j<X.columns(); ++j)
386 14 Aug 05 jari 60         X(i,j)=X(i,j-1)*x(i);
386 14 Aug 05 jari 61     md_.fit(X,y);
729 05 Jan 07 peter 62     chisq_ = md_.chisq();
383 12 Aug 05 jari 63   }
383 12 Aug 05 jari 64
713 21 Dec 06 peter 65
1120 21 Feb 08 peter 66   const utility::Vector& Polynomial::fit_parameters(void) const
713 21 Dec 06 peter 67   {
4200 19 Aug 22 peter 68     return md_.fit_parameters();
713 21 Dec 06 peter 69   }
713 21 Dec 06 peter 70
713 21 Dec 06 peter 71
586 19 Jun 06 peter 72   double Polynomial::predict(const double x) const
586 19 Jun 06 peter 73   {
1120 21 Feb 08 peter 74     utility::Vector vec(power_+1,1);
586 19 Jun 06 peter 75     for (size_t i=1; i<=power_; ++i)
586 19 Jun 06 peter 76       vec(i) = vec(i-1)*x;
586 19 Jun 06 peter 77     return md_.predict(vec);
586 19 Jun 06 peter 78   }
383 12 Aug 05 jari 79
586 19 Jun 06 peter 80
728 04 Jan 07 peter 81   double Polynomial::s2(void) const
728 04 Jan 07 peter 82   {
728 04 Jan 07 peter 83     return chisq()/(ap_.n()-power_-1);
728 04 Jan 07 peter 84   }
728 04 Jan 07 peter 85
728 04 Jan 07 peter 86
727 04 Jan 07 peter 87   double Polynomial::standard_error2(const double x) const
586 19 Jun 06 peter 88   {
1120 21 Feb 08 peter 89     utility::Vector vec(power_+1,1);
586 19 Jun 06 peter 90     for (size_t i=1; i<=power_; ++i)
586 19 Jun 06 peter 91       vec(i) = vec(i-1)*x;
727 04 Jan 07 peter 92     return md_.standard_error2(vec);
586 19 Jun 06 peter 93   }
586 19 Jun 06 peter 94
681 11 Oct 06 jari 95 }}} // of namespaces regression, yat, and theplu