source: trunk/yat/regression/LinearWeighted.cc

Last change on this file was 2919, checked in by Peter, 9 years ago

update copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.0 KB
Line 
1// $Id: LinearWeighted.cc 2919 2012-12-19 06:54:23Z peter $
2
3/*
4  Copyright (C) 2005 Peter Johansson
5  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
6  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2012 Peter Johansson
8
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include <config.h>
26
27#include "LinearWeighted.h"
28#include "yat/statistics/AveragerPairWeighted.h"
29#include "yat/utility/Vector.h"
30
31#include <cassert>
32
33namespace theplu {
34namespace yat {
35namespace regression {
36
37  LinearWeighted::LinearWeighted(void)
38    : OneDimensionalWeighted(), alpha_(0), alpha_var_(0), beta_(0),
39      beta_var_(0)
40  {
41  }
42
43  LinearWeighted::~LinearWeighted(void)
44  {
45  }
46
47  double LinearWeighted::alpha(void) const
48  {
49    return alpha_;
50  }
51
52  double LinearWeighted::alpha_var(void) const
53  {
54    return alpha_var_;
55  }
56
57  double LinearWeighted::beta(void) const
58  {
59    return beta_;
60  }
61
62  double LinearWeighted::beta_var(void) const
63  {
64    return beta_var_;
65  }
66
67  void LinearWeighted::fit(const utility::VectorBase& x,
68                           const utility::VectorBase& y,
69                           const utility::VectorBase& w)
70  {
71    assert(x.size()==y.size());
72    assert(x.size()==w.size());
73
74    // AveragerPairWeighted requires 2 weights but works only on the
75    // product wx*wy, so we can send in w and a dummie to get what we
76    // want.
77    ap_.reset();
78    yat::utility::Vector dummy(x.size(), 1.0);
79    add(ap_, x.begin(), x.end(), y.begin(),dummy.begin(),w.begin());
80
81    alpha_ = m_y();
82    beta_ = sxy()/sxx();
83
84    chisq_=0;
85    for (size_t i=0; i<x.size(); ++i){
86      double res = predict(x(i))-y(i);
87      chisq_ += w(i)*res*res;
88    }
89
90    alpha_var_ = s2()/ap_.y_averager().sum_w();
91    beta_var_ = s2()/sxx(); 
92  }
93
94  double LinearWeighted::m_x(void) const
95  {
96    return ap_.x_averager().mean();
97  }
98
99  double LinearWeighted::m_y(void) const
100  {
101    return ap_.y_averager().mean();
102  }
103
104  double LinearWeighted::predict(const double x) const
105  { 
106    return alpha_ + beta_ * (x-m_x()); 
107  }
108
109 
110  double LinearWeighted::s2(double w) const
111  {
112    return chisq_/(w*(ap_.y_averager().n()-2));
113  }
114
115  double LinearWeighted::standard_error2(const double x) const
116  {
117    return alpha_var_ + beta_var_*(x-m_x())*(x-m_x());
118  }
119
120
121  double LinearWeighted::sxx(void) const
122  {
123    return ap_.x_averager().sum_xx_centered();
124  }
125
126
127  double LinearWeighted::sxy(void) const
128  {
129    return ap_.sum_xy_centered();
130  }
131
132
133  double LinearWeighted::syy(void) const
134  {
135    return ap_.y_averager().sum_xx_centered();
136  }
137
138}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.