Changeset 216


Ignore:
Timestamp:
Dec 29, 2004, 10:29:54 AM (18 years ago)
Author:
Peter
Message:

RegressionLocal? added

Location:
trunk/src
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Makefile.am

    r207 r216  
    1414  Kernel.cc kNNI.cc matrix.cc Merge.cc NNI.cc PCA.cc \
    1515  PolynomialKernelFunction.cc random_singleton.cc \
     16  RegressionKernelBox.cc \
    1617  RegressionLinear.cc RegressionLocal.cc RegressionNaive.cc ROC.cc \
    1718  Score.cc Statistics.cc \
    … …  
    3233  KernelFunction.h kNNI.h matrix.h Merge.h NNI.h PCA.h Pearson.h \
    3334  PolynomialKernelFunction.h random_singleton.h \
    34   Regression.h RegressionLinear.h RegressionLocal.h RegressionNaive.h \
    35   ROC.h Score.h \
     35  Regression.h RegressionKernel.h RegressionKernelBox.h \
     36  RegressionLinear.h RegressionLocal.h RegressionNaive.h ROC.h Score.h \
    3637  PolynomialKernelFunction.h random_singleton.h Regression.h \
    37   RegressionKernel.h RegressionLinear.h ROC.h Score.h \
     38  RegressionKernel.h RegressionKernelBox.h \
     39  RegressionLinear.h ROC.h Score.h \
    3840  Statistics.h stl_utility.h SVD.h SVM.h tScore.h vector.h WeNNI.h \
    3941  WeightedAverager.h

  • trunk/src/Regression.h

    r213 r216  
    3434    virtual ~Regression(void) {};
    3535         
     36    virtual void estimate(const double x, double& y, double& y_err) const=0;
     37
    3638    ///
    3739    /// This function computes the best-fit given a model (see
    … …  
    4446    /// This function computes the best-fit given a model (see
    4547    /// specific class for details) by minimizing \f$
    46     /// \sum{(w_i\hat{y_i}-y_i)^2} \f$, where \f$ \hat{y} \f$ is the
     48    /// \sum{w_i(\hat{y_i}-y_i)^2} \f$, where \f$ \hat{y} \f$ is the
    4749    /// fitted value. The weight \f$ w_i \f$ is should be proportional
    4850    /// to the inverse of the variance for \f$ y_i \f$

  • trunk/src/RegressionKernel.h

    r205 r216  
    88
    99  ///
    10   ///
     10  /// Virtual class for calculating the weights in Local Regression
    1111  ///
    1212
    1313  class RegressionKernel
    1414  {
    15     RegressionKernel(void);
     15   
     16    public:
     17    ///
     18    /// Constructor
     19    ///
     20    RegressionKernel();
     21
     22    ///
     23    /// Function calculating the weight
     24    ///
     25    virtual double weight(const double) const=0;
    1626  };
    1727

  • trunk/src/RegressionLinear.cc

    r213 r216  
    2020  }
    2121
    22 
    2322}} // of namespace statistics and namespace theplu

  • trunk/src/RegressionLinear.h

    r213 r216  
    4040    virtual ~RegressionLinear(void) {};
    4141         
     42    inline void estimate(const double x, double& y, double& y_err) const
     43    {
     44      gsl_fit_linear_est(x, m_, k_, cov00_, cov01_, cov11_, &y, &y_err);
     45    }
     46
    4247    ///
    4348    /// This function computes the best-fit linear regression
    … …  
    4752    ///
    4853    inline int fit(const gslapi::vector& x, const gslapi::vector& y)
    49     { return gsl_fit_linear(x.gsl_vector_pointer()->data,
     54    {
     55      return gsl_fit_linear(x.gsl_vector_pointer()->data,
    5056                            x.gsl_vector_pointer()->stride,
    5157                            y.gsl_vector_pointer()->data,
    … …  
    6773    inline int fit_weighted(const gslapi::vector& x, const gslapi::vector& y,
    6874                            const gslapi::vector& w)
    69     { return gsl_fit_wlinear(x.gsl_vector_pointer()->data,
     75    {
     76      return gsl_fit_wlinear(x.gsl_vector_pointer()->data,
    7077                             x.gsl_vector_pointer()->stride,
    7178                             w.gsl_vector_pointer()->data,

  • trunk/src/RegressionLocal.cc

    r206 r216  
    88#include "vector.h"
    99
     10#include "algorithm"
    1011
    1112namespace theplu {
    … …  
    1415  RegressionLocal::RegressionLocal(const gslapi::vector& x,
    1516                                   const gslapi::vector& y,
    16                                    const Regression& r,
    17                                    const RegressionKernel& k)
    18     : kernel_(&k), regression_(&r)
     17                                   Regression& r,
     18                                   RegressionKernel& k,
     19                                   const gslapi::vector& estimated_x)
     20    : kernel_(&k), regression_(&r), estimated_x_(estimated_x),
     21      estimated_y_(estimated_x), estimated_y_err_(estimated_x)
    1922  {
     23    //sorting data with respect to x
     24    // Peter, could this be done without all the copying
     25    std::vector<std::pair<double, double> > data;
     26    for (size_t i=0; i<x.size(); i++){
     27      std::pair<double, double> tmp(x(i),y(i));
     28      data_.push_back(tmp);
     29    }
     30    sort(data_.begin(), data_.end());
     31  }
     32
     33  void RegressionLocal::fit(const double f)
     34  {
     35    // number of points on each side
     36    size_t nof_data_points = static_cast<size_t>(f/2*data_.size());
     37    size_t index_min = 0;
     38    size_t index_max = 0;
     39    for (size_t i=0; i<estimated_x_.size(); i++) {
     40      // extreme left case
     41      if (data_[nof_data_points-1].first>estimated_x_(i)){
     42        index_min = 0;
     43        index_max = 2*nof_data_points-1;
     44      }
     45     
     46      // extreme right case
     47      else if (data_[data_.size()-estimated_x_.size()-1].first<
     48               estimated_x_(i)){
     49        index_min = data_.size()-2*nof_data_points;
     50        index_max = data_.size()-1;
     51      }
     52     
     53      // normal case - window in the middle
     54      else {
     55        for (size_t j=nof_data_points; j<data_.size(); j++)
     56          if (data_[j].first>estimated_x_(i)) {
     57            index_min = j-nof_data_points;
     58            index_max = j+nof_data_points-1;
     59            break;
     60          }
     61      }
     62     
     63      // copying data
     64      // Peter, this is a stupid solution which could be
     65      // improved when the vector class is updated.
     66      gslapi::vector x(index_max-index_min+1);
     67      gslapi::vector y(index_max-index_min+1);
     68      gslapi::vector w(index_max-index_min+1);
     69      for (size_t j=0; j<x.size(); j++)
     70        x(j)=data_[index_min+j].first;
     71      for (size_t j=0; j<y.size(); j++)
     72        y(j)=data_[index_min+j].second;
     73      // calculating weights
     74      for (size_t j=0; j<y.size(); j++)
     75        w(j) = kernel_->weight( 2*(x(j)-data_[index_min].first)/
     76                                (data_[index_max].first-
     77                                 data_[index_min].first)-1);
     78
     79     
     80      // fitting the regressor locally
     81      regression_->fit_weighted(x,y,w);
     82      regression_->estimate(estimated_x_(i),
     83                            estimated_y_(i),
     84                            estimated_y_err_(i));
     85    }
    2086  }
    2187
    2288
     89  std::ostream& operator<< (std::ostream& s, const RegressionLocal& r)
     90  {
     91    gslapi::vector x=r.estimated_x();
     92    gslapi::vector y(r.estimated_y());
     93    gslapi::vector y_err(r.estimated_y_err());
     94    for (size_t i=0; i<x.size(); i++){
     95      s << x(i) << "\t";
     96      s << y(i) << "\t";
     97      s << y_err(i) << "\n";
     98    }
     99    return s;
     100  }
     101
    23102}} // of namespace cpptools and namespace theplu

  • trunk/src/RegressionLocal.h

    r206 r216  
    3232
    3333    ///
    34     /// Constructor doing the thing.
     34    /// Constructor loading the object with data, type of regressor
     35    /// and type kernel.
    3536    ///
    36     RegressionLocal(const gslapi::vector& x, const gslapi::vector& y,
    37                     const Regression&, const RegressionKernel&);
     37    RegressionLocal(const gslapi::vector& x, const gslapi::vector& y,
     38                    Regression& r, RegressionKernel& k,
     39                    const gslapi::vector&);
    3840
    3941    ///
    … …  
    4749    virtual ~RegressionLocal(void) {};
    4850
     51    inline gslapi::vector estimated_x(void) const { return estimated_x_; }
     52    inline gslapi::vector estimated_y(void) const { return estimated_y_; }
     53    inline gslapi::vector estimated_y_err(void) const {return estimated_y_err_;}
     54 
     55    ///
     56    /// Function
     57    ///
     58    void fit(const double fraction);
    4959         
    50   private:
    51     gslapi::vector data_x_;
     60  private:
     61    std::vector<std::pair<double, double> > data_;
    5262    gslapi::vector data_y_;
    53     const RegressionKernel* kernel_;
    54     const Regression* regression_;
     63    RegressionKernel* kernel_;
     64    Regression* regression_;
    5565    gslapi::vector estimated_x_;
    5666    gslapi::vector estimated_y_;
     67    gslapi::vector estimated_y_err_;
    5768   
     69   
    5870  };
     71
     72  ///
     73  /// The output operator for Local Regression Class.
     74  ///
     75  std::ostream& operator<< (std::ostream& s, const RegressionLocal&);
    5976
    6077}} // of namespace statistics and namespace theplu

  • trunk/src/RegressionNaive.h

    r213 r216  
    1313////////////////////////
    1414#include <gsl/gsl_fit.h>
     15#include <utility>
    1516
    1617namespace theplu {
    … …  
    4041    ///
    4142    virtual ~RegressionNaive(void) {};
     43
     44    inline void estimate(const double x, double& y, double& y_err)
     45    { y=m_; y_err=sqrt(var_); }
    4246         
    4347    ///
Note: See TracChangeset for help on using the changeset viewer.