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

// $Id: Local.cc 430 2005-12-08 22:53:08Z 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>
#include <iostream>

namespace theplu {
namespace statistics {
namespace regression {


  void Local::fit(const double width, const gslapi::vector& points)
  {
    // Peter, avoid this copying
    x_=points;
    y_predicted_ = gslapi::vector(points.size());
    y_err_ = gslapi::vector(points.size());

    sort(data_.begin(), data_.end());

    // coying data to 2 gslapi vectors ONCE to use views from
    gslapi::vector x(data_.size());
    gslapi::vector y(data_.size());
    for (size_t j=0; j<x.size(); j++)
      x(j)=data_[j].first;
    for (size_t j=0; j<y.size(); j++)
      y(j)=data_[j].second;
    
    size_t min_index = 0;
    size_t max_index = 1; 

    for (size_t i=0; i<points.size(); i++) {
      
      while (points(min_index) < x(i)-width)
        min_index++;
      if (min_index)
        min_index--;
      while (points(max_index) < x(i)+width && max_index<points.size()-1)
        max_index++;

      gslapi::vector x_local(x, min_index, max_index-min_index+1);
      gslapi::vector y_local(y, min_index, max_index-min_index+1);

      // calculating weights
      gslapi::vector w(max_index-min_index+1);
      for (size_t j=0; j<y_local.size(); j++)
        w(j) = kernel_->weight( (x_local(j)- x(i))/width );
      
      // fitting the regressor locally
      regressor_->fit(x_local,y_local,w);
      regressor_->predict(x(i), y_predicted_(i),y_err_(i));
    }
  }

  std::ostream& operator<<(std::ostream& os, const Local& r)
  {
    os << "# column 1: x\n"
      << "# column 2: y\n"
      << "# column 3: y_err\n";
    for (size_t i=0; i<r.x().size(); i++) {
      os << r.x()(i) << "\t" 
         << r.y_predicted()(i) << "\t"
         << r.y_err()(i) << "\n";
    }   

    return os;
  }
}}} // of namespaces regression, statisitcs and thep