source: trunk/c++_tools/statistics/LinearWeighted.cc @ 586

Last change on this file since 586 was 586, checked in by Peter, 15 years ago

closes #23 redesign of regression classes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 1.2 KB
Line 
1// $Id: LinearWeighted.cc 586 2006-06-19 09:56:04Z peter $
2
3#include <c++_tools/statistics/LinearWeighted.h>
4
5#include <c++_tools/statistics/AveragerPairWeighted.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 LinearWeighted::fit(const gslapi::vector& x,
17                           const gslapi::vector& y,
18                           const gslapi::vector& w)
19  {
20    // AveragerPairWeighted requires 2 weights but works only on the
21    // product wx*wy, so we can send in w and a dummie to get what we
22    // want.
23    gslapi::vector dummie(w.size(),1);
24    AveragerPairWeighted ap;
25    ap.add_values(x,y,w,dummie);
26
27    double m_x = ap.x_averager().mean();
28    double m_y = ap.y_averager().mean();
29   
30    double sxy = ap.sum_xy_centered();
31
32    double sxx = ap.x_averager().sum_xx_centered();
33    double syy = ap.y_averager().sum_xx_centered();
34
35    // estimating the noise level. see attached document for motivation
36    // of the expression.
37    s2_= (syy-sxy*sxy/sxx)/(w.sum()-2*(w*w)/w.sum()) ;
38   
39    alpha_ = m_y;
40    beta_ = sxy/sxx;
41    alpha_var_ = ap.y_averager().standard_error() * 
42      ap.y_averager().standard_error();
43    beta_var_ = s2_/sxx; 
44    m_x_=m_x;
45  }
46
47
48}}} // of namespaces regression, statisitcs and thep
Note: See TracBrowser for help on using the repository browser.