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

Last change on this file since 429 was 429, checked in by Peter, 18 years ago

separating weighted and non-weighted regression to different classes.

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 2.4 KB
Line
1// \$Id: Local.cc 429 2005-12-08 19:50:11Z 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
11namespace theplu {
12namespace statistics {
13namespace regression {
14
15
16  void Local::fit(std::ostream& s, const double f,
17                            const u_int step_size)
18  {
19    sort(data_.begin(), data_.end());
20
21    // coying data to 2 gslapi vectors ONCE to use views from
22    gslapi::vector x_total(data_.size());
23    gslapi::vector y_total(data_.size());
24    for (size_t j=0; j<x_total.size(); j++)
25      x_total(j)=data_[j].first;
26    for (size_t j=0; j<y_total.size(); j++)
27      y_total(j)=data_[j].second;
28
30    s << std::endl;
31
32    // number of points on each side
33    size_t nof_points = static_cast<size_t>(f*data_.size());
34    double width = 0;
35    int min_index = 0;
36    for (size_t i=0; i<data_.size(); i+=step_size) {
37      // determining boundaries for window/kernel
38
39      // extreme left case
40      if (data_[i].first <
41          (data_[nof_points-1].first+data_[nof_points-1].first)/2 ){
42        min_index=0;
43        width = data_[nof_points-1].first-data_[0].first;
44      }
45      // extreme right case
46      else if (data_[i].first > (data_[data_.size()-1].first+
47                            data_[data_.size()-nof_points].first)/2 ){
48        min_index = data_.size() - nof_points;
49        width = data_[i].first - data_[data_.size()-nof_points].first;
50      }
51      else {
52        // stepping forward until x is in the middle
53        while(data_[min_index+nof_points-1].first-data_[i].first
54              < data_[i].first-data_[min_index].first)
55          min_index++;
56        width = data_[min_index+nof_points-1].first-data_[min_index].first;
57        // Peter, should check if x gets closer to centrum if we step
58        // back one step
59      }
60
61
62      gslapi::vector x(x_total, min_index, nof_points);
63      gslapi::vector y(y_total, min_index, nof_points);
64
65      // calculating weights
66      gslapi::vector w(nof_points);
67      for (size_t j=0; j<y.size(); j++)
68        w(j) = kernel_->weight( (x(j)- (data_[min_index].first +
69                                        data_[min_index+nof_points-1].first) /2)
70                                /width );
71
72      // fitting the regressor locally
73      regressor_->fit(x,y,w);
74      x_.push_back(data_[i].first);
75      y_.push_back(0.0);
76      y_err_.push_back(0.0);
77      regressor_->predict(x_[x_.size()-1], y_[y_.size()-1],
78                          y_err_[y_err_.size()-1]);
79      regressor_->print(s);
80      s << std::endl;
81    }
82  }
83
84
85}}} // of namespaces regression, statisitcs and thep
Note: See TracBrowser for help on using the repository browser.