source: trunk/src/RegressionLocal.cc @ 235

// $Id: RegressionLocal.cc 235 2005-02-21 14:53:48Z peter $

#include "RegressionLocal.h"

#include "Regression.h"
#include "RegressionKernel.h"
#include "RegressionLinear.h"
#include "vector.h"

#include "algorithm"

namespace theplu {
namespace statistics {

  RegressionLocal::RegressionLocal(Regression& r, 
                                   RegressionKernel& k)
    : kernel_(&k), regressor_(&r)
  {
  }

  void RegressionLocal::fit(const double f, const u_int step_size)
  {
    sort(data_.begin(), data_.end());

    // number of points on each side
    size_t nof_points = static_cast<size_t>(f*data_.size());
    double width = 0;
    int min_index = 0;
    for (size_t i=0; i<data_.size(); i+=step_size) {
      x_.push_back(data_[i].first);
      
      // determining boundaries window/kernel
      
      // extreme left case
      if (data_[i].first <
          (data_[nof_points-1].first+data_[nof_points-1].first)/2 ){
        min_index=0;
        width = data_[nof_points-1].first-data_[i].first;
      }
      // extreme right case
      else if (data_[i].first > (data_[data_.size()-1].first+
                            data_[data_.size()-nof_points].first)/2 ){
        min_index = data_.size() - nof_points;
        width = data_[i].first - data_[data_.size()-nof_points].first;
      }
      else {
        // stepping forward until x is in the middle
        while(data_[min_index+nof_points-1].first-data_[i].first 
              < data_[i].first-data_[min_index].first)
          min_index++;
        width = data_[min_index+nof_points-1].first-data_[i].first;
        // Peter, should check if x gets closer to centrum if we step
        // back one step
      }


      // copying data 
      // Peter, too much copying. Move the copying outside loop and
      // use a vector view (when available in gslapi::vector class).
      gslapi::vector x(nof_points);
      gslapi::vector y(nof_points);
      gslapi::vector w(nof_points);
      for (size_t j=0; j<x.size(); j++)
        x(j)=data_[min_index+j].first;
      for (size_t j=0; j<y.size(); j++){
        y(j)=data_[min_index+j].second;
      }
      // calculating weights
      for (size_t j=0; j<y.size(); j++)
        w(j) = kernel_->weight( (x(j)-(data_[min_index].first+
                                       data_[min_index+nof_points-1].first) /2)
                                /width );
      
      // fitting the regressor locally
      regressor_->fit(x,y,w);
      regressor_->predict(x_[i], y_[i], y_err_[i]);
    }
  }

  std::ostream& RegressionLocal::print(std::ostream& s) const
  {
    for (size_t i=0; i<x_.size(); i++) {
      regressor_->print(s);
      s << std::endl;
    }
    return s;
  }


}} // of namespace cpptools and namespace theplu