Changeset 235 for trunk/src/RegressionLocal.cc

Timestamp:

Feb 21, 2005, 3:53:48 PM (18 years ago)

Author:

Peter

Message:

Major modifications

File:

• trunk/src/RegressionLocal.cc (modified) (1 diff)

trunk/src/RegressionLocal.cc

@@ -13 +13 @@
 namespace statistics {
 
-  RegressionLocal::RegressionLocal(const gslapi::vector& x, 
-                                   const gslapi::vector& y, 
-                                   Regression& r, 
-                                   RegressionKernel& k,
-                                   const size_t nof_predictions)
-    : kernel_(&k), regression_(&r)
+  RegressionLocal::RegressionLocal(Regression& r, 
+                                   RegressionKernel& k)
+    : kernel_(&k), regressor_(&r)
   {
-    estimated_x_ = gslapi::vector(nof_predictions);
-    estimated_y_ = estimated_x_;
-    estimated_y_err_ = estimated_x_;
-    //sorting data with respect to x
-    // Peter, could this be done without all the copying
-    std::vector<std::pair<double, double> > data;
-    for (size_t i=0; i<x.size(); i++){
-      std::pair<double, double> tmp(x(i),y(i));
-      data_.push_back(tmp);
-    }
-    sort(data_.begin(), data_.end());
-    for (size_t i=0; i<nof_predictions; i++)
-      estimated_x_(i) =
-        data_[static_cast<size_t>((i+0.5)/nof_predictions*data_.size())].first;
   }
 
-  void RegressionLocal::fit(const double f)
+  void RegressionLocal::fit(const double f, const u_int step_size)
   {
+    sort(data_.begin(), data_.end());
+
     // number of points on each side
-    size_t nof_data_points = static_cast<size_t>(f/2*data_.size());
-    size_t index_min = 0;
-    size_t index_max = 0;
-    for (size_t i=0; i<estimated_x_.size(); i++) {
+    size_t nof_points = static_cast<size_t>(f*data_.size());
+    double width = 0;
+    int min_index = 0;
+    for (size_t i=0; i<data_.size(); i+=step_size) {
+      x_.push_back(data_[i].first);
+      
+      // determining boundaries window/kernel
+      
       // extreme left case
-      if (data_[nof_data_points-1].first>estimated_x_(i)){
-        index_min = 0;
-        index_max = 2*nof_data_points-1;
+      if (data_[i].first <
+          (data_[nof_points-1].first+data_[nof_points-1].first)/2 ){
+        min_index=0;
+        width = data_[nof_points-1].first-data_[i].first;
       }
-      
       // extreme right case
-      else if (data_[data_.size()-estimated_x_.size()-1].first<
-               estimated_x_(i)){
-        index_min = data_.size()-2*nof_data_points;
-        index_max = data_.size()-1;
+      else if (data_[i].first > (data_[data_.size()-1].first+
+                            data_[data_.size()-nof_points].first)/2 ){
+        min_index = data_.size() - nof_points;
+        width = data_[i].first - data_[data_.size()-nof_points].first;
       }
-      
-      // normal case - window in the middle
       else {
-        for (size_t j=nof_data_points; j<data_.size(); j++)
-          if (data_[j].first>estimated_x_(i)) {
-            index_min = j-nof_data_points;
-            index_max = j+nof_data_points-1;
-            break;
-          }
+        // stepping forward until x is in the middle
+        while(data_[min_index+nof_points-1].first-data_[i].first 
+              < data_[i].first-data_[min_index].first)
+          min_index++;
+        width = data_[min_index+nof_points-1].first-data_[i].first;
+        // Peter, should check if x gets closer to centrum if we step
+        // back one step
       }
-      
+
+
       // copying data 
-      // Peter, this is a stupid solution which could be
-      // improved by using a moving view (in the vector class).
-      gslapi::vector x(index_max-index_min+1);
-      gslapi::vector y(index_max-index_min+1);
-      gslapi::vector w(index_max-index_min+1);
+      // Peter, too much copying. Move the copying outside loop and
+      // use a vector view (when available in gslapi::vector class).
+      gslapi::vector x(nof_points);
+      gslapi::vector y(nof_points);
+      gslapi::vector w(nof_points);
       for (size_t j=0; j<x.size(); j++)
-        x(j)=data_[index_min+j].first;
+        x(j)=data_[min_index+j].first;
       for (size_t j=0; j<y.size(); j++){
-        y(j)=data_[index_min+j].second;
-        if (y(j)>10)
-          std::cout << x(j) << "\t" << y(j) << std::endl;
+        y(j)=data_[min_index+j].second;
       }
       // calculating weights
       for (size_t j=0; j<y.size(); j++)
-        w(j) = kernel_->weight( 2*(x(j)-data_[index_min].first)/
-                                (data_[index_max].first-
-                                 data_[index_min].first)-1);
-
+        w(j) = kernel_->weight( (x(j)-(data_[min_index].first+
+                                       data_[min_index+nof_points-1].first) /2)
+                                /width );
       
       // fitting the regressor locally
-      std::cout << "x: " << estimated_x_(i) << std::endl;
-      std::cout << "x.begin: " << x(0) << std::endl;
-      std::cout << "x.end: " << x(x.size()-1) << std::endl;
-      regression_->fit(x,y,w);
-      std::cout << "estimating value\n";
-      regression_->predict(estimated_x_(i),
-                            estimated_y_(i),
-                            estimated_y_err_(i));
+      regressor_->fit(x,y,w);
+      regressor_->predict(x_[i], y_[i], y_err_[i]);
     }
   }
-  
-  std::ostream& operator<< (std::ostream& s, const RegressionLocal& r)
+
+  std::ostream& RegressionLocal::print(std::ostream& s) const
   {
-    gslapi::vector x=r.estimated_x();
-    gslapi::vector y(r.estimated_y());
-    gslapi::vector y_err(r.estimated_y_err());
-    for (size_t i=0; i<x.size(); i++){
-      s << x(i) << "\t";
-      s << y(i) << "\t";
-      s << y_err(i) << "\n";
+    for (size_t i=0; i<x_.size(); i++) {
+      regressor_->print(s);
+      s << std::endl;
     }
     return s;
   }
 
+
 }} // of namespace cpptools and namespace theplu