source: trunk/lib/statistics/Local.cc @ 429

// $Id: Local.cc 429 2005-12-08 19:50:11Z peter $

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

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

//#include <algorithm>

namespace theplu {
namespace statistics {
namespace regression {


  void Local::fit(std::ostream& s, const double f, 
                            const u_int step_size)
  {
    sort(data_.begin(), data_.end());

    // coying data to 2 gslapi vectors ONCE to use views from
    gslapi::vector x_total(data_.size());
    gslapi::vector y_total(data_.size());
    for (size_t j=0; j<x_total.size(); j++)
      x_total(j)=data_[j].first;
    for (size_t j=0; j<y_total.size(); j++)
      y_total(j)=data_[j].second;
    
    regressor_->print_header(s); 
    s << std::endl;

    // 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) {
      // determining boundaries for 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_[0].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_[min_index].first;
        // Peter, should check if x gets closer to centrum if we step
        // back one step
      }

      
      gslapi::vector x(x_total, min_index, nof_points);
      gslapi::vector y(y_total, min_index, nof_points);

      // calculating weights
      gslapi::vector w(nof_points);
      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);
      x_.push_back(data_[i].first);
      y_.push_back(0.0);
      y_err_.push_back(0.0);
      regressor_->predict(x_[x_.size()-1], y_[y_.size()-1], 
                          y_err_[y_err_.size()-1]);
      regressor_->print(s);
      s << std::endl;
    }
  }


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