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

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

Continued fixing of regression stuff.

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