yat/regression/MultiDimensionalWeighted.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 /*
2119 12 Dec 09 peter 4   Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 5   Copyright (C) 2012, 2020, 2022, 2023 Peter Johansson
586 19 Jun 06 peter 6
1437 25 Aug 08 peter 7   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 8
675 10 Oct 06 jari 9   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 10   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 11   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 12   License, or (at your option) any later version.
675 10 Oct 06 jari 13
675 10 Oct 06 jari 14   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 15   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 17   General Public License for more details.
675 10 Oct 06 jari 18
675 10 Oct 06 jari 19   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 20   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 21 */
675 10 Oct 06 jari 22
2881 18 Nov 12 peter 23 #include <config.h>
2881 18 Nov 12 peter 24
680 11 Oct 06 jari 25 #include "MultiDimensionalWeighted.h"
739 12 Jan 07 peter 26 #include "yat/statistics/AveragerWeighted.h"
1121 22 Feb 08 peter 27 #include "yat/utility/Matrix.h"
1120 21 Feb 08 peter 28 #include "yat/utility/Vector.h"
675 10 Oct 06 jari 29
586 19 Jun 06 peter 30 #include <cassert>
586 19 Jun 06 peter 31
586 19 Jun 06 peter 32 namespace theplu {
680 11 Oct 06 jari 33 namespace yat {
586 19 Jun 06 peter 34 namespace regression {
586 19 Jun 06 peter 35
703 18 Dec 06 jari 36   MultiDimensionalWeighted::MultiDimensionalWeighted(void)
703 18 Dec 06 jari 37     : chisquare_(0), work_(NULL)
703 18 Dec 06 jari 38   {
703 18 Dec 06 jari 39   }
703 18 Dec 06 jari 40
703 18 Dec 06 jari 41   MultiDimensionalWeighted::~MultiDimensionalWeighted(void)
703 18 Dec 06 jari 42   {
3907 09 May 20 peter 43     gsl_multifit_linear_free(work_);
703 18 Dec 06 jari 44   }
703 18 Dec 06 jari 45
741 13 Jan 07 peter 46
741 13 Jan 07 peter 47   double MultiDimensionalWeighted::chisq() const
741 13 Jan 07 peter 48   {
741 13 Jan 07 peter 49     return chisquare_;
741 13 Jan 07 peter 50   }
741 13 Jan 07 peter 51
741 13 Jan 07 peter 52
4125 14 Jan 22 peter 53   void MultiDimensionalWeighted::fit(const utility::MatrixBase& x,
1022 01 Feb 08 peter 54                                      const utility::VectorBase& y,
1022 01 Feb 08 peter 55                                      const utility::VectorBase& w)
586 19 Jun 06 peter 56   {
586 19 Jun 06 peter 57     assert(y.size()==w.size());
586 19 Jun 06 peter 58     assert(x.rows()==y.size());
586 19 Jun 06 peter 59
1098 18 Feb 08 jari 60     covariance_.resize(x.columns(),x.columns());
3907 09 May 20 peter 61     fit_parameters_.resize(x.columns());
3907 09 May 20 peter 62     gsl_multifit_linear_free(work_);
750 17 Feb 07 jari 63     if (!(work_=gsl_multifit_linear_alloc(x.rows(),fit_parameters_.size())))
750 17 Feb 07 jari 64       throw utility::GSL_error("MultiDimensionalWeighted::fit failed to allocate memory");
750 17 Feb 07 jari 65     int status = gsl_multifit_wlinear(x.gsl_matrix_p(), w.gsl_vector_p(),
4200 19 Aug 22 peter 66                                       y.gsl_vector_p(),
750 17 Feb 07 jari 67                                       fit_parameters_.gsl_vector_p(),
750 17 Feb 07 jari 68                                       covariance_.gsl_matrix_p(), &chisquare_,
750 17 Feb 07 jari 69                                       work_);
750 17 Feb 07 jari 70     if (status)
750 17 Feb 07 jari 71       throw utility::GSL_error(std::string("MultiDimensionalWeighted::fit",
750 17 Feb 07 jari 72                                            status));
739 12 Jan 07 peter 73
739 12 Jan 07 peter 74     statistics::AveragerWeighted aw;
916 30 Sep 07 peter 75     add(aw, y.begin(), y.end(), w.begin());
739 12 Jan 07 peter 76     s2_ = chisquare_ / (aw.n()-fit_parameters_.size());
739 12 Jan 07 peter 77     covariance_ *= s2_;
586 19 Jun 06 peter 78   }
586 19 Jun 06 peter 79
740 12 Jan 07 peter 80
1120 21 Feb 08 peter 81   const utility::Vector& MultiDimensionalWeighted::fit_parameters(void) const
740 12 Jan 07 peter 82   {
740 12 Jan 07 peter 83     return fit_parameters_;
740 12 Jan 07 peter 84   }
740 12 Jan 07 peter 85
789 10 Mar 07 jari 86
1022 01 Feb 08 peter 87   double MultiDimensionalWeighted::predict(const utility::VectorBase& x) const
718 26 Dec 06 jari 88   {
731 06 Jan 07 peter 89     assert(x.size()==fit_parameters_.size());
718 26 Dec 06 jari 90     return fit_parameters_ * x;
718 26 Dec 06 jari 91   }
718 26 Dec 06 jari 92
789 10 Mar 07 jari 93
1022 01 Feb 08 peter 94   double MultiDimensionalWeighted::prediction_error2(const utility::VectorBase& x,
739 12 Jan 07 peter 95                                                      const double w) const
586 19 Jun 06 peter 96   {
740 12 Jan 07 peter 97     return standard_error2(x) + s2(w);
586 19 Jun 06 peter 98   }
586 19 Jun 06 peter 99
586 19 Jun 06 peter 100
740 12 Jan 07 peter 101   double MultiDimensionalWeighted::s2(const double w) const
740 12 Jan 07 peter 102   {
740 12 Jan 07 peter 103     return s2_/w;
740 12 Jan 07 peter 104   }
740 12 Jan 07 peter 105
740 12 Jan 07 peter 106
4200 19 Aug 22 peter 107   double
1022 01 Feb 08 peter 108   MultiDimensionalWeighted::standard_error2(const utility::VectorBase& x) const
586 19 Jun 06 peter 109   {
4292 04 Feb 23 peter 110     return x * covariance_ * x;
586 19 Jun 06 peter 111   }
586 19 Jun 06 peter 112
681 11 Oct 06 jari 113 }}} // of namespaces regression, yat, and theplu

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