source: trunk/src/RegressionNaive.h @ 221

// $Id: RegressionNaive.h 221 2004-12-30 22:36:25Z peter $

#ifndef _theplu_statistics_regression_naive_
#define _theplu_statistics_regression_naive_

// C++ tools include
/////////////////////
#include "Averager.h"
#include "Regression.h"
#include "vector.h"
#include "WeightedAverager.h"
// Standard C++ includes
////////////////////////
//#include <gsl/gsl_fit.h>
#include <utility>

namespace theplu {
namespace statistics {  
  
  ///
  /// Class for Regression.   
  ///
  
  class RegressionNaive : public Regression 
  {
  
  public:
    ///
    /// Default Constructor.
    ///
    RegressionNaive(void);

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

    ///
    /// Destructor 
    ///
    virtual ~RegressionNaive(void) {};
         
    ///
    /// This function computes the best-fit for the naive model \f$
    /// y = m \f$ from vectors \a x and \a y, by minimizing \f$
    /// \sum{(y_i-m)^2} \f$.
    ///
    inline void fit(const gslapi::vector& x, const gslapi::vector& y)
    {
      Averager a;
      for (size_t i=0; i<y.size(); i++)
        a.add(y(i));
      m_=a.mean();
      s2_=a.variance();
      m_err_=a.standard_error();
    }

    ///
    /// This function computes the best-fit for the naive model \f$ y
    /// = m \f$ from vectors \a x and \a y, by minimizing \f$ \sum
    /// w_i(y_i-m)^2 \f$. The weight \f$ w_i \f$ is proportional to
    /// the inverse of the variance for \f$ y_i \f$
    ///
    inline void fit(const gslapi::vector& x,
                    const gslapi::vector& y,
                    const gslapi::vector& w)
    {
      WeightedAverager a;
      for (size_t i=0; i<y.size(); i++)
        a.add(y(i), w(i));
      m_=a.mean();
      m_err_=a.standard_error();
      s2_=m_err_*m_err_*w.sum(); 
    }

    ///
    /// @return the parameter m
    ///
    inline double m(void) const { return m_; }

    ///
    /// @return standard deviation of parameter m
    ///
    inline double alpha_var(void) const { return m_err_; }

    ///
    /// Function predicting value using the naive model. \a y_err is
    /// the expected deviation from the line for a new data point. If
    /// weights are used a weight can be specified for the new point.
    ///
    inline void predict(const double x, double& y, double& y_err, 
                        const double w=1.0) const 
    { 
      y = m_; 
      y_err = sqrt(m_err_*m_err_ + s2_/w);
    }

         
  private:
    double s2_; // noise level ( var = s2/w in weighted case)
    double m_;
    double m_err_;
  };

}} // of namespace statistics and namespace theplu

#endif