source: trunk/src/RegressionLinear.cc @ 221

// $Id: RegressionLinear.cc 221 2004-12-30 22:36:25Z peter $

// Header
#include "RegressionLinear.h"

// C++_tools
#include "AveragerPair.h"
#include "Regression.h"
#include "vector.h"

#include <gsl/gsl_fit.h>


namespace theplu {
namespace statistics {

  
  RegressionLinear::RegressionLinear(void)
    : Regression(), alpha_(0), alpha_var_(0), beta_(0), beta_var_(0), m_x_(0),
      s2_(0)
  {
  }

  void RegressionLinear::fit(const gslapi::vector& x, const gslapi::vector& y) 
  {
    statistics::AveragerPair ap;
    for (size_t i=0; i<x.size(); i++)
      ap.add(x(i),y(i));

    alpha_ = ap.y_averager().mean();
    beta_ = ap.covariance() / ap.x_averager().variance();

    // estimating the noise level, i.e. the conditional variance of y
    // given x, Var(y|x).
    double Q = (ap.y_averager().sum_xsqr_centered() - ap.sum_xy_centered() *
                ap.sum_xy_centered()/ap.x_averager().sum_xsqr_centered() ) ;
    s2_ = Q/(x.size()-2);
    r2_= 1;
    alpha_var_ = s2_ / x.size();
    beta_var_ = s2_ / ap.x_averager().sum_xsqr_centered();
    m_x_ = ap.x_averager().mean();
  }

  void RegressionLinear::fit(const gslapi::vector& x, const gslapi::vector& y,
                             const gslapi::vector& w)
  {
    double m_x = w*x /w.sum();
    double m_y = w*y /w.sum();
    
    double sxy = 0;
    for (size_t i=0; i<x.size(); i++) 
      sxy += w(i)*(x(i)-m_x)*(y(i)-m_y);

    double sxx = 0;
    for (size_t i=0; i<x.size(); i++) 
      sxx += w(i)*(x(i)-m_x)*(x(i)-m_x);

    double syy = 0;
    for (size_t i=0; i<y.size(); i++) 
      syy += w(i)*(y(i)-m_y)*(y(i)-m_y);

    // estimating the noise level. see attached document for motivation
    // of the expression.
    s2_= (syy-sxy*sxy/sxx)/(w.sum()-2*(w*w)/w.sum()) ;
    
    alpha_ = m_y;
    beta_ = sxy/sxx;
    alpha_var_ = s2_/w.sum();
    beta_var_ = s2_/sxx;  
    m_x_=m_x;
  }

    
}} // of namespace statistics and namespace theplu