source: trunk/c++_tools/statistics/Local.cc @ 586

Last change on this file since 586 was 586, checked in by Peter, 15 years ago

closes #23 redesign of regression classes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.9 KB
Line 
1// $Id: Local.cc 586 2006-06-19 09:56:04Z peter $
2
3#include <c++_tools/statistics/Local.h>
4
5#include <c++_tools/gslapi/vector.h>
6#include <c++_tools/statistics/Kernel.h>
7#include <c++_tools/statistics/OneDimensionalWeighted.h>
8
9#include <algorithm>
10#include <cassert>
11#include <iostream>
12
13namespace theplu {
14namespace statistics {
15namespace regression {
16
17
18  void Local::fit(const size_t step_size, const size_t nof_points)
19  {
20    if (step_size==0 || nof_points<3){
21      // Peter to Jari, throw exception?
22      std::cerr << "theplu::statistics::regression::Local "
23                << "Parameters invalid. Fitting ignored." << std::endl;
24      return;
25    }
26
27    size_t nof_fits=data_.size()/step_size;
28    x_= gslapi::vector(nof_fits);
29    y_predicted_ = gslapi::vector(x_.size());
30    y_err_ = gslapi::vector(x_.size());
31    sort(data_.begin(), data_.end());
32
33    // coying data to 2 gslapi vectors ONCE to use views from
34    gslapi::vector x(data_.size());
35    gslapi::vector y(data_.size());
36    for (size_t j=0; j<x.size(); j++){
37      x(j)=data_[j].first;
38      y(j)=data_[j].second;
39    }
40
41    // looping over regression points and perform local regression
42    for (size_t i=0; i<nof_fits; i++) {
43      size_t max_index = static_cast<size_t>( (i+0.5)*step_size );
44      size_t min_index;
45      double width; // distance from middle of windo to border of window
46      double x_mid; // middle of window
47      // right border case
48      if (max_index > data_.size()-1){
49        min_index = max_index - nof_points + 1;
50        max_index = data_.size()-1;
51        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
52                  ( 2*(max_index-min_index)) );
53        x_mid = x(min_index)+width;
54      }
55      // normal middle case
56      else if (max_index > nof_points-1){
57        min_index = max_index - nof_points + 1;
58        width = (x(max_index)-x(min_index))/2;
59        x_mid = x(min_index)+width;
60      }
61      // left border case
62      else {
63        min_index = 0;
64        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
65                  ( 2*(max_index-min_index)) );
66        x_mid = x(max_index)-width;
67      }
68      assert(min_index<data_.size());
69      assert(max_index<data_.size());
70                               
71      gslapi::vector x_local(x, min_index, max_index-min_index+1);
72      gslapi::vector y_local(y, min_index, max_index-min_index+1);
73
74      // calculating weights
75      gslapi::vector w(max_index-min_index+1);
76      for (size_t j=0; j<w.size(); j++)
77        w(j) = kernel_->weight( (x_local(j)- x_mid)/width );
78     
79      // fitting the regressor locally
80      regressor_->fit(x_local,y_local,w);
81      assert(i<y_predicted_.size());
82      assert(i<y_err_.size());
83      y_predicted_(i) = regressor_->predict(x(i*step_size));
84      y_err_(i) = regressor_->standard_error(x(i*step_size));
85    }
86  }
87
88  std::ostream& operator<<(std::ostream& os, const Local& r)
89  {
90    os << "# column 1: x\n"
91      << "# column 2: y\n"
92      << "# column 3: y_err\n";
93    for (size_t i=0; i<r.x().size(); i++) {
94      os << r.x()(i) << "\t" 
95         << r.y_predicted()(i) << "\t"
96         << r.y_err()(i) << "\n";
97    }   
98
99    return os;
100  }
101}}} // of namespaces regression, statisitcs and thep
Note: See TracBrowser for help on using the repository browser.