source:trunk/lib/statistics/RegressionLocal.cc@303

Last change on this file since 303 was 303, checked in by Peter, 17 years ago

docs

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