source: trunk/lib/statistics/Linear.cc @ 383

Last change on this file since 383 was 383, checked in by Jari Häkkinen, 16 years ago

Moved regression stuff into a namespace ...statistics::regression.
Filenames was changed and a first draft of a Polynomial regression is
implemented (compiles but does not run properly).

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