source: trunk/src/RegressionLinear.cc @ 235

Last change on this file since 235 was 235, checked in by Peter, 17 years ago

Major modifications

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.4 KB
Line 
1// $Id: RegressionLinear.cc 235 2005-02-21 14:53:48Z peter $
2
3// Header
4#include "RegressionLinear.h"
5
6// C++_tools
7#include "AveragerPair.h"
8#include "Regression.h"
9#include "vector.h"
10
11#include <gsl/gsl_fit.h>
12
13
14namespace theplu {
15namespace statistics {
16
17 
18  RegressionLinear::RegressionLinear(void)
19    : Regression(), alpha_(0), alpha_var_(0), beta_(0), beta_var_(0), m_x_(0),
20      s2_(0)
21  {
22  }
23
24  void RegressionLinear::fit(const gslapi::vector& x, const gslapi::vector& y) 
25  {
26    statistics::AveragerPair ap;
27    for (size_t i=0; i<x.size(); i++)
28      ap.add(x(i),y(i));
29
30    alpha_ = ap.y_averager().mean();
31    beta_ = ap.covariance() / ap.x_averager().variance();
32
33    // estimating the noise level, i.e. the conditional variance of y
34    // given x, Var(y|x).
35    double Q = (ap.y_averager().sum_xsqr_centered() - ap.sum_xy_centered() *
36                ap.sum_xy_centered()/ap.x_averager().sum_xsqr_centered() ) ;
37    s2_ = Q/(x.size()-2);
38    r2_= 1;
39    alpha_var_ = s2_ / x.size();
40    beta_var_ = s2_ / ap.x_averager().sum_xsqr_centered();
41    m_x_ = ap.x_averager().mean();
42  }
43
44  void RegressionLinear::fit(const gslapi::vector& x, const gslapi::vector& y,
45                             const gslapi::vector& w)
46  {
47    double m_x = w*x /w.sum();
48    double m_y = w*y /w.sum();
49   
50    double sxy = 0;
51    for (size_t i=0; i<x.size(); i++) 
52      sxy += w(i)*(x(i)-m_x)*(y(i)-m_y);
53
54    double sxx = 0;
55    for (size_t i=0; i<x.size(); i++) 
56      sxx += w(i)*(x(i)-m_x)*(x(i)-m_x);
57
58    double syy = 0;
59    for (size_t i=0; i<y.size(); i++) 
60      syy += w(i)*(y(i)-m_y)*(y(i)-m_y);
61
62    // estimating the noise level. see attached document for motivation
63    // of the expression.
64    s2_= (syy-sxy*sxy/sxx)/(w.sum()-2*(w*w)/w.sum()) ;
65   
66    alpha_ = m_y;
67    beta_ = sxy/sxx;
68    alpha_var_ = s2_/w.sum();
69    beta_var_ = s2_/sxx; 
70    m_x_=m_x;
71  }
72
73  std::ostream& RegressionLinear::print(std::ostream& s) const
74  {
75    s << x_ << "\t"
76      << y_ << "\t"
77      << y_err_ << "\t"
78      << alpha_ << "\t"
79      << alpha_err() << "\t"
80      << beta_ << "\t"
81      << beta_err() << "\t"
82      << s2_ << "\t"
83      << r2_;
84   
85    return s;
86  }
87   
88  std::ostream& RegressionLinear::print_header(std::ostream& s) const
89  {
90    s << "# column 1: x\n"
91      << "# column 2: y\n"
92      << "# column 3: y_err\n"
93      << "# column 4: alpha\n"
94      << "# column 5: alpha_err\n"
95      << "# column 6: beta\n"
96      << "# column 7: beta_err\n"
97      << "# column 8: s2 (var(y|x))\n"
98      << "# column 9: r2 (coefficient of determination)";
99    return s;
100  }
101
102   
103}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.