source: trunk/src/RegressionLocal.cc @ 216

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

RegressionLocal? added

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.8 KB
Line 
1// $Id: RegressionLocal.cc 216 2004-12-29 09:29:54Z peter $
2
3#include "RegressionLocal.h"
4
5#include "Regression.h"
6#include "RegressionKernel.h"
7#include "RegressionLinear.h"
8#include "vector.h"
9
10#include "algorithm"
11
12namespace theplu {
13namespace statistics {
14
15  RegressionLocal::RegressionLocal(const gslapi::vector& x, 
16                                   const gslapi::vector& y, 
17                                   Regression& r, 
18                                   RegressionKernel& k,
19                                   const gslapi::vector& estimated_x)
20    : kernel_(&k), regression_(&r), estimated_x_(estimated_x), 
21      estimated_y_(estimated_x), estimated_y_err_(estimated_x)
22  {
23    //sorting data with respect to x
24    // Peter, could this be done without all the copying
25    std::vector<std::pair<double, double> > data;
26    for (size_t i=0; i<x.size(); i++){
27      std::pair<double, double> tmp(x(i),y(i));
28      data_.push_back(tmp);
29    }
30    sort(data_.begin(), data_.end());
31  }
32
33  void RegressionLocal::fit(const double f)
34  {
35    // number of points on each side
36    size_t nof_data_points = static_cast<size_t>(f/2*data_.size());
37    size_t index_min = 0;
38    size_t index_max = 0;
39    for (size_t i=0; i<estimated_x_.size(); i++) {
40      // extreme left case
41      if (data_[nof_data_points-1].first>estimated_x_(i)){
42        index_min = 0;
43        index_max = 2*nof_data_points-1;
44      }
45     
46      // extreme right case
47      else if (data_[data_.size()-estimated_x_.size()-1].first<
48               estimated_x_(i)){
49        index_min = data_.size()-2*nof_data_points;
50        index_max = data_.size()-1;
51      }
52     
53      // normal case - window in the middle
54      else {
55        for (size_t j=nof_data_points; j<data_.size(); j++)
56          if (data_[j].first>estimated_x_(i)) {
57            index_min = j-nof_data_points;
58            index_max = j+nof_data_points-1;
59            break;
60          }
61      }
62     
63      // copying data
64      // Peter, this is a stupid solution which could be
65      // improved when the vector class is updated.
66      gslapi::vector x(index_max-index_min+1);
67      gslapi::vector y(index_max-index_min+1);
68      gslapi::vector w(index_max-index_min+1);
69      for (size_t j=0; j<x.size(); j++)
70        x(j)=data_[index_min+j].first;
71      for (size_t j=0; j<y.size(); j++)
72        y(j)=data_[index_min+j].second;
73      // calculating weights
74      for (size_t j=0; j<y.size(); j++)
75        w(j) = kernel_->weight( 2*(x(j)-data_[index_min].first)/
76                                (data_[index_max].first-
77                                 data_[index_min].first)-1);
78
79     
80      // fitting the regressor locally
81      regression_->fit_weighted(x,y,w);
82      regression_->estimate(estimated_x_(i),
83                            estimated_y_(i),
84                            estimated_y_err_(i));
85    }
86  }
87
88
89  std::ostream& operator<< (std::ostream& s, const RegressionLocal& r)
90  {
91    gslapi::vector x=r.estimated_x();
92    gslapi::vector y(r.estimated_y());
93    gslapi::vector y_err(r.estimated_y_err());
94    for (size_t i=0; i<x.size(); i++){
95      s << x(i) << "\t";
96      s << y(i) << "\t";
97      s << y_err(i) << "\n";
98    }
99    return s;
100  }
101
102}} // of namespace cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.