Changeset 430


Ignore:
Timestamp:
Dec 8, 2005, 11:53:08 PM (18 years ago)
Author:
Peter
Message:

changed interface of Regression::Local

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/lib/statistics/Linear.cc

    r429 r430  
    3434  }
    3535
    36   void Linear::predict(const double x, double& y, double& y_err)
     36  void Linear::predict(const double x, double& y, double& y_err) const
    3737  {
    3838    y = alpha_ + beta_ * (x-m_x_);

  • trunk/lib/statistics/Linear.h

    r429 r430  
    6767    /// the expected deviation from the line for a new data point.
    6868    ///
    69     void predict(const double x, double& y, double& y_err);
     69    void predict(const double x, double& y, double& y_err) const;
    7070
    7171    ///

  • trunk/lib/statistics/LinearWeighted.h

    r429 r430  
    7575    {
    7676      y = alpha_ + beta_ * (x-m_x_);
    77       y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)+s2_/w );
     77      y_err_ = y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)+s2_/w );
    7878      x_=x;
    7979      y_=y;
    80       y_err_=y_err;
    8180    }
    8281

  • trunk/lib/statistics/Local.cc

    r429 r430  
    77#include <c++_tools/statistics/OneDimensionalWeighted.h>
    88
    9 //#include <algorithm>
     9#include <algorithm>
     10#include <iostream>
    1011
    1112namespace theplu {
    … …  
    1415
    1516
    16   void Local::fit(std::ostream& s, const double f,
    17                             const u_int step_size)
     17  void Local::fit(const double width, const gslapi::vector& points)
    1818  {
     19    // Peter, avoid this copying
     20    x_=points;
     21    y_predicted_ = gslapi::vector(points.size());
     22    y_err_ = gslapi::vector(points.size());
     23
    1924    sort(data_.begin(), data_.end());
    2025
    2126    // coying data to 2 gslapi vectors ONCE to use views from
    22     gslapi::vector x_total(data_.size());
    23     gslapi::vector y_total(data_.size());
    24     for (size_t j=0; j<x_total.size(); j++)
    25       x_total(j)=data_[j].first;
    26     for (size_t j=0; j<y_total.size(); j++)
    27       y_total(j)=data_[j].second;
     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;
    2833   
    29     regressor_->print_header(s);
    30     s << std::endl;
     34    size_t min_index = 0;
     35    size_t max_index = 1;
    3136
    32     // number of points on each side
    33     size_t nof_points = static_cast<size_t>(f*data_.size());
    34     double width = 0;
    35     int min_index = 0;
    36     for (size_t i=0; i<data_.size(); i+=step_size) {
    37       // determining boundaries for window/kernel
     37    for (size_t i=0; i<points.size(); i++) {
    3838     
    39       // extreme left case
    40       if (data_[i].first <
    41           (data_[nof_points-1].first+data_[nof_points-1].first)/2 ){
    42         min_index=0;
    43         width = data_[nof_points-1].first-data_[0].first;
    44       }
    45       // extreme right case
    46       else if (data_[i].first > (data_[data_.size()-1].first+
    47                             data_[data_.size()-nof_points].first)/2 ){
    48         min_index = data_.size() - nof_points;
    49         width = data_[i].first - data_[data_.size()-nof_points].first;
    50       }
    51       else {
    52         // stepping forward until x is in the middle
    53         while(data_[min_index+nof_points-1].first-data_[i].first
    54               < data_[i].first-data_[min_index].first)
    55           min_index++;
    56         width = data_[min_index+nof_points-1].first-data_[min_index].first;
    57         // Peter, should check if x gets closer to centrum if we step
    58         // back one step
    59       }
     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++;
    6045
    61      
    62       gslapi::vector x(x_total, min_index, nof_points);
    63       gslapi::vector y(y_total, min_index, nof_points);
     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);
    6448
    6549      // calculating weights
    66       gslapi::vector w(nof_points);
    67       for (size_t j=0; j<y.size(); j++)
    68         w(j) = kernel_->weight( (x(j)- (data_[min_index].first +
    69                                         data_[min_index+nof_points-1].first) /2)
    70                                 /width );
     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 );
    7153     
    7254      // fitting the regressor locally
    73       regressor_->fit(x,y,w);
    74       x_.push_back(data_[i].first);
    75       y_.push_back(0.0);
    76       y_err_.push_back(0.0);
    77       regressor_->predict(x_[x_.size()-1], y_[y_.size()-1],
    78                           y_err_[y_err_.size()-1]);
    79       regressor_->print(s);
    80       s << std::endl;
     55      regressor_->fit(x_local,y_local,w);
     56      regressor_->predict(x(i), y_predicted_(i),y_err_(i));
    8157    }
    8258  }
    8359
     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    }   
    8470
     71    return os;
     72  }
    8573}}} // of namespaces regression, statisitcs and thep

  • trunk/lib/statistics/Local.h

    r429 r430  
    88#include <c++_tools/gslapi/vector.h>
    99
     10#include <iostream>
    1011
    1112namespace theplu {
    … …  
    4849
    4950    ///
    50     /// Function returning the points where to predict
    51     ///
    52     inline const std::vector<double>& x(void) const { return x_; }
    53 
    54     ///
    5551    /// Function returning predicted values
    5652    ///
    57     inline const std::vector<double>& y(void) const { return y_; }
     53    inline const gslapi::vector& y_predicted(void) const
     54    { return y_predicted_; }
    5855
    5956    ///
    6057    /// Function returning error of predictions
    6158    ///
    62     inline const std::vector<double>& y_err(void) const { return y_err_; }
     59    inline const gslapi::vector& y_err(void) const { return y_err_; }
    6360 
    6461    ///
    6562    /// Performs the fit in data defined by add using a
    6663    /// RegressionKernel and a Regression method defined in the
    67     /// constructor. The function starts by regressing over the lowest
    68     /// x value, followed by stepping forward \a step_size number of
    69     /// points, et cetera. The kernel is centralized over the x-value
    70     /// in question and the width is set so \a fraction of all points
    71     /// are used. The result is sent to ostream \a s, using print
    72     /// function in used Regression.
     64    /// constructor. For each element in vector \a x a fit is
     65    /// performed. The kernel used has width \a width and is
     66    /// centralized over the data point \f$ x_i \f$, which means data
     67    /// in \f$ [x-width, x+width] is used in the fit.
    7368    ///
    74     void fit(std::ostream& s, const double fraction, const u_int step_size=1);
     69    void fit(const double width, const gslapi::vector& x);
    7570
     71    ///
     72    /// @return x-values where fitting was performed.
     73    ///
     74    inline const gslapi::vector& x(void) const { return x_; }
    7675
    7776  private:
    … …  
    8483    Kernel* kernel_;
    8584    OneDimensionalWeighted* regressor_;
    86     std::vector<double> x_;
    87     std::vector<double> y_;
    88     std::vector<double> y_err_;
     85    gslapi::vector x_;
     86    gslapi::vector y_predicted_;
     87    gslapi::vector y_err_;
    8988
    9089  };
    9190
    92 }}} // of namespaces regression, statisitcs and thep
     91  ///
     92  /// The output operator for the RegressionLocal class.
     93  ///
     94  std::ostream& operator<<(std::ostream&, const Local& );
     95
     96
     97}}} // of namespaces regression, statistics and thep
    9398
    9499#endif

  • trunk/lib/statistics/Naive.cc

    r429 r430  
    2323  }
    2424
    25   void Naive::predict(const double x, double& y, double& y_err)
     25  void Naive::predict(const double x, double& y, double& y_err) const
    2626  {
    2727    y = ap_.x_averager().mean();

  • trunk/lib/statistics/Naive.h

    r429 r430  
    5454    /// estimation of the mean.
    5555    ///
    56     void predict(const double x, double& y, double& y_err) ;
     56    void predict(const double x, double& y, double& y_err) const;
    5757
    5858    ///

  • trunk/test/regression_test.cc

    r429 r430  
    4242  gslapi::vector w(4); w(0)=0.1;  w(1)=0.2;  w(2)=0.3;  w(3)=0.4;
    4343
    44   // testing regression::LinearWeighted
     44  *error << "testing regression::LinearWeighted" << std::endl;
    4545  statistics::regression::LinearWeighted linear_w;
    4646  linear_w.fit(x,y,w);
    … …  
    4949  linear_w.predict(1990,y_predicted,y_predicted_err);
    5050  if (y_predicted!=12.8){
    51     *error << "regression_Linear: cannot reproduce fit." << std::endl;
     51    *error << "regression_LinearWeighted: cannot reproduce fit." << std::endl;
    5252    ok=false;
    5353  }
    5454
    5555  // testing regression::NaiveWeighted
     56  *error << "testing regression::NaiveWeighted" << std::endl;
    5657  statistics::regression::NaiveWeighted naive_w;
    5758  naive_w.fit(x,y,w);
    … …  
    6566
    6667  // testing regression::Local
     68  *error << "testing regression::Local" << std::endl;
    6769  statistics::regression::KernelBox kb;
    6870  statistics::regression::LinearWeighted rl;
    … …  
    7880
    7981  // testing regression::Polynomial
     82  *error << "testing regression::Polynomial" << std::endl;
    8083  {
    8184    std::ifstream s("data/regression_gauss.data");
    … …  
    98101  }
    99102
     103  *error << "testing regression::Linear" << std::endl;
     104  statistics::regression::Linear lin;
     105 
     106  *error << "testing regression::Naive" << std::endl;
     107  statistics::regression::Naive naive;
     108
    100109  if (error!=&std::cerr)
    101110    delete error;
    … …  
    110119{
    111120  statistics::regression::Local rl(r,k);
     121  gslapi::vector points(1000);
    112122  for (size_t i=0; i<1000; i++){
    113     double x = i;
    114     double y = 10.0;
    115     rl.add(x, y);
     123    rl.add(i, 10);
     124    points(i)=i;
    116125  }
    117126
    118   std::ofstream myout("/dev/null");
    119   rl.fit(myout,0.1,1);
     127  rl.fit(100, points);
    120128
    121   std::vector<double> y = rl.y();
     129  gslapi::vector y = rl.y_predicted();
    122130  for (size_t i=0; i<y.size(); i++)
    123     if (y[i]!=10.0){
     131    if (y(i)!=10.0){
    124132      return false;
    125133    }
Note: See TracChangeset for help on using the changeset viewer.