source: trunk/yat/regression/LinearWeighted.cc @ 1703

Last change on this file since 1703 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.0 KB
Line 
1// $Id: LinearWeighted.cc 1487 2008-09-10 08:41:36Z jari $
2
3/*
4  Copyright (C) 2005 Peter Johansson
5  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
6  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2008 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 "LinearWeighted.h"
26#include "yat/statistics/AveragerPairWeighted.h"
27#include "yat/utility/Vector.h"
28
29#include <cassert>
30
31namespace theplu {
32namespace yat {
33namespace regression {
34
35  LinearWeighted::LinearWeighted(void)
36    : OneDimensionalWeighted(), alpha_(0), alpha_var_(0), beta_(0),
37      beta_var_(0)
38  {
39  }
40
41  LinearWeighted::~LinearWeighted(void)
42  {
43  }
44
45  double LinearWeighted::alpha(void) const
46  {
47    return alpha_;
48  }
49
50  double LinearWeighted::alpha_var(void) const
51  {
52    return alpha_var_;
53  }
54
55  double LinearWeighted::beta(void) const
56  {
57    return beta_;
58  }
59
60  double LinearWeighted::beta_var(void) const
61  {
62    return beta_var_;
63  }
64
65  void LinearWeighted::fit(const utility::VectorBase& x,
66                           const utility::VectorBase& y,
67                           const utility::VectorBase& w)
68  {
69    assert(x.size()==y.size());
70    assert(x.size()==w.size());
71
72    // AveragerPairWeighted requires 2 weights but works only on the
73    // product wx*wy, so we can send in w and a dummie to get what we
74    // want.
75    ap_.reset();
76    yat::utility::Vector dummy(x.size(), 1.0);
77    add(ap_, x.begin(), x.end(), y.begin(),dummy.begin(),w.begin());
78
79    alpha_ = m_y();
80    beta_ = sxy()/sxx();
81
82    chisq_=0;
83    for (size_t i=0; i<x.size(); ++i){
84      double res = predict(x(i))-y(i);
85      chisq_ += w(i)*res*res;
86    }
87
88    alpha_var_ = s2()/ap_.y_averager().sum_w();
89    beta_var_ = s2()/sxx(); 
90  }
91
92  double LinearWeighted::m_x(void) const
93  {
94    return ap_.x_averager().mean();
95  }
96
97  double LinearWeighted::m_y(void) const
98  {
99    return ap_.y_averager().mean();
100  }
101
102  double LinearWeighted::predict(const double x) const
103  { 
104    return alpha_ + beta_ * (x-m_x()); 
105  }
106
107 
108  double LinearWeighted::s2(double w) const
109  {
110    return chisq_/(w*(ap_.y_averager().n()-2));
111  }
112
113  double LinearWeighted::standard_error2(const double x) const
114  {
115    return alpha_var_ + beta_var_*(x-m_x())*(x-m_x());
116  }
117
118
119  double LinearWeighted::sxx(void) const
120  {
121    return ap_.x_averager().sum_xx_centered();
122  }
123
124
125  double LinearWeighted::sxy(void) const
126  {
127    return ap_.sum_xy_centered();
128  }
129
130
131  double LinearWeighted::syy(void) const
132  {
133    return ap_.y_averager().sum_xx_centered();
134  }
135
136}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.