Changeset 235 for trunk/src/RegressionLocal.cc
 Timestamp:
 Feb 21, 2005, 3:53:48 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/RegressionLocal.cc
r221 r235 13 13 namespace statistics { 14 14 15 RegressionLocal::RegressionLocal(const gslapi::vector& x, 16 const gslapi::vector& y, 17 Regression& r, 18 RegressionKernel& k, 19 const size_t nof_predictions) 20 : kernel_(&k), regression_(&r) 15 RegressionLocal::RegressionLocal(Regression& r, 16 RegressionKernel& k) 17 : kernel_(&k), regressor_(&r) 21 18 { 22 estimated_x_ = gslapi::vector(nof_predictions);23 estimated_y_ = estimated_x_;24 estimated_y_err_ = estimated_x_;25 //sorting data with respect to x26 // Peter, could this be done without all the copying27 std::vector<std::pair<double, double> > data;28 for (size_t i=0; i<x.size(); i++){29 std::pair<double, double> tmp(x(i),y(i));30 data_.push_back(tmp);31 }32 sort(data_.begin(), data_.end());33 for (size_t i=0; i<nof_predictions; i++)34 estimated_x_(i) =35 data_[static_cast<size_t>((i+0.5)/nof_predictions*data_.size())].first;36 19 } 37 20 38 void RegressionLocal::fit(const double f )21 void RegressionLocal::fit(const double f, const u_int step_size) 39 22 { 23 sort(data_.begin(), data_.end()); 24 40 25 // number of points on each side 41 size_t nof_data_points = static_cast<size_t>(f/2*data_.size()); 42 size_t index_min = 0; 43 size_t index_max = 0; 44 for (size_t i=0; i<estimated_x_.size(); i++) { 26 size_t nof_points = static_cast<size_t>(f*data_.size()); 27 double width = 0; 28 int min_index = 0; 29 for (size_t i=0; i<data_.size(); i+=step_size) { 30 x_.push_back(data_[i].first); 31 32 // determining boundaries window/kernel 33 45 34 // extreme left case 46 if (data_[nof_data_points1].first>estimated_x_(i)){ 47 index_min = 0; 48 index_max = 2*nof_data_points1; 35 if (data_[i].first < 36 (data_[nof_points1].first+data_[nof_points1].first)/2 ){ 37 min_index=0; 38 width = data_[nof_points1].firstdata_[i].first; 49 39 } 50 51 40 // extreme right case 52 else if (data_[ data_.size()estimated_x_.size()1].first<53 estimated_x_(i)){54 index_min = data_.size()2*nof_data_points;55 index_max = data_.size()1;41 else if (data_[i].first > (data_[data_.size()1].first+ 42 data_[data_.size()nof_points].first)/2 ){ 43 min_index = data_.size()  nof_points; 44 width = data_[i].first  data_[data_.size()nof_points].first; 56 45 } 57 58 // normal case  window in the middle59 46 else { 60 for (size_t j=nof_data_points; j<data_.size(); j++) 61 if (data_[j].first>estimated_x_(i)) { 62 index_min = jnof_data_points; 63 index_max = j+nof_data_points1; 64 break; 65 } 47 // stepping forward until x is in the middle 48 while(data_[min_index+nof_points1].firstdata_[i].first 49 < data_[i].firstdata_[min_index].first) 50 min_index++; 51 width = data_[min_index+nof_points1].firstdata_[i].first; 52 // Peter, should check if x gets closer to centrum if we step 53 // back one step 66 54 } 67 55 56 68 57 // copying data 69 // Peter, t his is a stupid solution which could be70 // improved by using a moving view (in thevector class).71 gslapi::vector x( index_maxindex_min+1);72 gslapi::vector y( index_maxindex_min+1);73 gslapi::vector w( index_maxindex_min+1);58 // Peter, too much copying. Move the copying outside loop and 59 // use a vector view (when available in gslapi::vector class). 60 gslapi::vector x(nof_points); 61 gslapi::vector y(nof_points); 62 gslapi::vector w(nof_points); 74 63 for (size_t j=0; j<x.size(); j++) 75 x(j)=data_[ index_min+j].first;64 x(j)=data_[min_index+j].first; 76 65 for (size_t j=0; j<y.size(); j++){ 77 y(j)=data_[index_min+j].second; 78 if (y(j)>10) 79 std::cout << x(j) << "\t" << y(j) << std::endl; 66 y(j)=data_[min_index+j].second; 80 67 } 81 68 // calculating weights 82 69 for (size_t j=0; j<y.size(); j++) 83 w(j) = kernel_>weight( 2*(x(j)data_[index_min].first)/ 84 (data_[index_max].first 85 data_[index_min].first)1); 86 70 w(j) = kernel_>weight( (x(j)(data_[min_index].first+ 71 data_[min_index+nof_points1].first) /2) 72 /width ); 87 73 88 74 // fitting the regressor locally 89 std::cout << "x: " << estimated_x_(i) << std::endl; 90 std::cout << "x.begin: " << x(0) << std::endl; 91 std::cout << "x.end: " << x(x.size()1) << std::endl; 92 regression_>fit(x,y,w); 93 std::cout << "estimating value\n"; 94 regression_>predict(estimated_x_(i), 95 estimated_y_(i), 96 estimated_y_err_(i)); 75 regressor_>fit(x,y,w); 76 regressor_>predict(x_[i], y_[i], y_err_[i]); 97 77 } 98 78 } 99 100 std::ostream& operator<< (std::ostream& s, const RegressionLocal& r)79 80 std::ostream& RegressionLocal::print(std::ostream& s) const 101 81 { 102 gslapi::vector x=r.estimated_x(); 103 gslapi::vector y(r.estimated_y()); 104 gslapi::vector y_err(r.estimated_y_err()); 105 for (size_t i=0; i<x.size(); i++){ 106 s << x(i) << "\t"; 107 s << y(i) << "\t"; 108 s << y_err(i) << "\n"; 82 for (size_t i=0; i<x_.size(); i++) { 83 regressor_>print(s); 84 s << std::endl; 109 85 } 110 86 return s; 111 87 } 112 88 89 113 90 }} // of namespace cpptools and namespace theplu
Note: See TracChangeset
for help on using the changeset viewer.