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

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

changed interface of Regression::Local

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 1.8 KB
Line
1// \$Id: Local.cc 430 2005-12-08 22:53:08Z 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 <iostream>
11
12namespace theplu {
13namespace statistics {
14namespace regression {
15
16
17  void Local::fit(const double width, const gslapi::vector& points)
18  {
19    // Peter, avoid this copying
20    x_=points;
21    y_predicted_ = gslapi::vector(points.size());
22    y_err_ = gslapi::vector(points.size());
23
24    sort(data_.begin(), data_.end());
25
26    // coying data to 2 gslapi vectors ONCE to use views from
27    gslapi::vector x(data_.size());
28    gslapi::vector y(data_.size());
29    for (size_t j=0; j<x.size(); j++)
30      x(j)=data_[j].first;
31    for (size_t j=0; j<y.size(); j++)
32      y(j)=data_[j].second;
33
34    size_t min_index = 0;
35    size_t max_index = 1;
36
37    for (size_t i=0; i<points.size(); i++) {
38
39      while (points(min_index) < x(i)-width)
40        min_index++;
41      if (min_index)
42        min_index--;
43      while (points(max_index) < x(i)+width && max_index<points.size()-1)
44        max_index++;
45
46      gslapi::vector x_local(x, min_index, max_index-min_index+1);
47      gslapi::vector y_local(y, min_index, max_index-min_index+1);
48
49      // calculating weights
50      gslapi::vector w(max_index-min_index+1);
51      for (size_t j=0; j<y_local.size(); j++)
52        w(j) = kernel_->weight( (x_local(j)- x(i))/width );
53
54      // fitting the regressor locally
55      regressor_->fit(x_local,y_local,w);
56      regressor_->predict(x(i), y_predicted_(i),y_err_(i));
57    }
58  }
59
60  std::ostream& operator<<(std::ostream& os, const Local& r)
61  {
62    os << "# column 1: x\n"
63      << "# column 2: y\n"
64      << "# column 3: y_err\n";
65    for (size_t i=0; i<r.x().size(); i++) {
66      os << r.x()(i) << "\t"
67         << r.y_predicted()(i) << "\t"
68         << r.y_err()(i) << "\n";
69    }
70
71    return os;
72  }
73}}} // of namespaces regression, statisitcs and thep
Note: See TracBrowser for help on using the repository browser.