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

// $Id: Linear.cc 430 2005-12-08 22:53:08Z peter $

#include <c++_tools/statistics/Linear.h>

#include <c++_tools/statistics/AveragerPair.h>
#include <c++_tools/gslapi/vector.h>

#include <gsl/gsl_fit.h>


namespace theplu {
namespace statistics {
namespace regression {


  void Linear::fit(const gslapi::vector& x, const gslapi::vector& y) 
  {
    ap_.reset();
    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 Linear::predict(const double x, double& y, double& y_err) const
  { 
    y = alpha_ + beta_ * (x-m_x_); 
    y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)+s2_ ); 
  }

  std::ostream& Linear::print_header(std::ostream& s) const
  {
    s << "# column 1: x\n"
      << "# column 2: y\n"
      << "# column 3: y_err\n"
      << "# column 4: alpha\n"
      << "# column 5: alpha_err\n"
      << "# column 6: beta\n"
      << "# column 7: beta_err\n"
      << "# column 8: s2 (var(y|x))\n"
      << "# column 9: r2 (coefficient of determination)";
    return s;
  }


}}} // of namespaces regression, statisitcs and thep