source: trunk/src/RegressionLinear.h @ 216

// $Id: RegressionLinear.h 216 2004-12-29 09:29:54Z peter $

#ifndef _theplu_statistics_regression_linear_
#define _theplu_statistics_regression_linear_

// C++ tools include
/////////////////////
#include "Regression.h"
#include "vector.h"

// Standard C++ includes
////////////////////////
#include <gsl/gsl_fit.h>

namespace theplu {
namespace statistics {  
  
  ///
  /// Class for Regression.   
  ///

  
  class RegressionLinear : public Regression 
  {
  
  public:
    ///
    /// Default Constructor.
    ///
    RegressionLinear(void);

    ///
    /// Copy Constructor. (not implemented)
    ///
    RegressionLinear(const RegressionLinear&);

    ///
    /// Destructor 
    ///
    virtual ~RegressionLinear(void) {};
         
    inline void estimate(const double x, double& y, double& y_err) const
    { 
      gsl_fit_linear_est(x, m_, k_, cov00_, cov01_, cov11_, &y, &y_err); 
    }

    ///
    /// This function computes the best-fit linear regression
    /// coefficients (k,m) of the model  \f$ y = kx + m \f$ from
    /// vectors \a x and \a y, by minimizing \f$ \sum{(kx_i+m-y_i)^2}
    /// \f$.
    ///
    inline int fit(const gslapi::vector& x, const gslapi::vector& y) 
    { 
      return gsl_fit_linear(x.gsl_vector_pointer()->data, 
                            x.gsl_vector_pointer()->stride,
                            y.gsl_vector_pointer()->data, 
                            y.gsl_vector_pointer()->stride,
                            x.gsl_vector_pointer()->size,
                            &m_, &k_, &cov00_, &cov01_, &cov11_, &sumsq_); }

    inline double k(void) const { return k_; }
    inline double m(void) const { return m_; }


    ///
    /// This function computes the best-fit linear regression
    /// coefficients (k,m) of the model \f$ y = kx + m \f$ from
    /// vectors \a x and \a y, by minimizing \f$ \sum w_i(kx_i+m-y_i)^2
    /// \f$. The weight \f$ w_i \f$ is the inverse of the variance for
    /// \f$ y_i \f$
    ///
    inline int fit_weighted(const gslapi::vector& x, const gslapi::vector& y,
                            const gslapi::vector& w)
    {
      return gsl_fit_wlinear(x.gsl_vector_pointer()->data, 
                             x.gsl_vector_pointer()->stride,
                             w.gsl_vector_pointer()->data, 
                             w.gsl_vector_pointer()->stride,
                             y.gsl_vector_pointer()->data, 
                             y.gsl_vector_pointer()->stride,
                             x.gsl_vector_pointer()->size,
                             &m_, &k_,  &cov00_, &cov01_, &cov11_, &sumsq_); }

         
  private:
    double cov00_;
    double cov01_;
    double cov11_;    
    double k_;
    double m_;

  };

}} // of namespace statistics and namespace theplu

#endif