Changeset 495


Ignore:
Timestamp:
Jan 11, 2006, 6:27:22 PM (16 years ago)
Author:
Peter
Message:

modified interface for Regression::Local

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/lib/statistics/LinearWeighted.cc

    r489 r495  
    3939    alpha_ = m_y;
    4040    beta_ = sxy/sxx;
    41     alpha_var_ = s2_/w.sum();
     41    alpha_var_ = ap.y_averager().standard_error() *
     42      ap.y_averager().standard_error();
    4243    beta_var_ = s2_/sxx; 
    4344    m_x_=m_x;
  • trunk/lib/statistics/LinearWeighted.h

    r430 r495  
    6868    ///
    6969    /// Function predicting value using the linear model. \a y_err is
    70     /// the expected deviation from the line for a new data point. If
    71     /// weights are used a weight can be specified for the new point.
     70    /// the estimated error of y(x)
    7271    ///
    7372    inline void predict(const double x, double& y, double& y_err,
     
    7574    {
    7675      y = alpha_ + beta_ * (x-m_x_);
    77       y_err_ = y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)+s2_/w );
     76      y_err_ = y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_));
    7877      x_=x;
    7978      y_=y;
  • trunk/lib/statistics/Local.cc

    r430 r495  
    1515
    1616
    17   void Local::fit(const double width, const gslapi::vector& points)
     17  void Local::fit(const size_t step_size, const size_t nof_points)
    1818  {
    19     // Peter, avoid this copying
    20     x_=points;
    21     y_predicted_ = gslapi::vector(points.size());
    22     y_err_ = gslapi::vector(points.size());
     19    if (step_size==0 || nof_points<3){
     20      // Peter to Jari, throw exception?
     21      std::cerr << "theplu::statistics::regression::Local "
     22                << "Parameters invalid. Fitting ignored." << std::endl;
     23      return;
     24    }
    2325
     26    size_t nof_fits=data_.size()/step_size;
     27    x_= gslapi::vector(nof_fits);
     28    y_predicted_ = gslapi::vector(x_.size());
     29    y_err_ = gslapi::vector(x_.size());
    2430    sort(data_.begin(), data_.end());
    2531
     
    2733    gslapi::vector x(data_.size());
    2834    gslapi::vector y(data_.size());
    29     for (size_t j=0; j<x.size(); j++)
     35    for (size_t j=0; j<x.size(); j++){
    3036      x(j)=data_[j].first;
    31     for (size_t j=0; j<y.size(); j++)
    3237      y(j)=data_[j].second;
    33    
    34     size_t min_index = 0;
    35     size_t max_index = 1;
     38    }
    3639
    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 
     40    // looping over regression points and perform local regression
     41    for (size_t i=0; i<nof_fits; i++) {
     42      size_t max_index = static_cast<size_t>( (i+0.5)*step_size );
     43      size_t min_index;
     44      double width; // distance from middle of windo to border of window
     45      double x_mid; // middle of window
     46      // right border case
     47      if (max_index > data_.size()-1){
     48        min_index = max_index - nof_points + 1;
     49        max_index = data_.size()-1;
     50        width = ( (( x(max_index)-x(0) )*(nof_points-1)) /
     51                  ( 2*(max_index-min_index)) );
     52        x_mid = x(min_index)+width;
     53      }
     54      // normal middle case
     55      else if (max_index > nof_points-1){
     56        min_index = max_index - nof_points + 1;
     57        width = (x(max_index)-x(min_index))/2;
     58        x_mid = x(min_index)+width;
     59      }
     60      // left border case
     61      else {
     62        min_index = 0;
     63        width = ( (( x(max_index)-x(0) )*(nof_points-1)) /
     64                  ( 2*(max_index-min_index)) );
     65        x_mid = x(max_index)-width;
     66      }
     67      assert(min_index<data_.size());
     68      assert(max_index<data_.size());
     69                               
    4670      gslapi::vector x_local(x, min_index, max_index-min_index+1);
    4771      gslapi::vector y_local(y, min_index, max_index-min_index+1);
     
    4973      // calculating weights
    5074      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 );
     75      for (size_t j=0; j<w.size(); j++)
     76        w(j) = kernel_->weight( (x_local(j)- x_mid)/width );
    5377     
    5478      // fitting the regressor locally
    5579      regressor_->fit(x_local,y_local,w);
    56       regressor_->predict(x(i), y_predicted_(i),y_err_(i));
     80      assert(i<y_predicted_.size());
     81      assert(i<y_err_.size());
     82      regressor_->predict(x(i*step_size), y_predicted_(i),y_err_(i));
    5783    }
    5884  }
  • trunk/lib/statistics/Local.h

    r430 r495  
    6060 
    6161    ///
    62     /// Performs the fit in data defined by add using a
    63     /// RegressionKernel and a Regression method defined in the
    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.
     62    /// @todo doc
     63    /// @param nof_points Number of points used in each fit
     64    /// @param step_size Size of step between each fit
    6865    ///
    69     void fit(const double width, const gslapi::vector& x);
     66    void fit(const size_t step_size, const size_t nof_points);
    7067
    7168    ///
  • trunk/lib/statistics/NaiveWeighted.cc

    r429 r495  
    3232    x_ = x;
    3333    y = m_;
    34     y_err = sqrt(m_err_*m_err_ + s2_/w);
     34    y_err = m_err_;
    3535  }
    3636 
  • trunk/lib/statistics/NaiveWeighted.h

    r429 r495  
    4848
    4949    ///
    50     /// Function predicting value using the naive model. \a y_err is
    51     /// the expected deviation from the line for a new data point. The
    52     /// weight for the new point can be specified. A smaller weight
    53     /// means larger error. The error has two components: the variance
    54     /// of point and error in estimation of m_.
     50    /// Function predicting value using the naive model, i.e. a
     51    /// weighted average. \a y_err is the standard error of the
     52    /// weighted mean
     53    ///
     54    /// @see AveragerWeighted::standard_error
    5555    ///
    5656    void predict(const double x, double& y, double& y_err,
  • trunk/lib/statistics/OneDimensionalWeighted.h

    r429 r495  
    4444                    const gslapi::vector& w)=0;
    4545    ///
    46     /// function predicting in one point
     46    /// function predicting in one point. y will contain the y(x)
     47    /// according to the model. y_err will contain the estimated error
     48    /// of y(x).
    4749    ///
    4850    virtual void predict(const double x, double& y, double& y_err,
  • trunk/test/regression_test.cc

    r489 r495  
    121121{
    122122  statistics::regression::Local rl(r,k);
    123   gslapi::vector points(1000);
    124123  for (size_t i=0; i<1000; i++){
    125124    rl.add(i, 10);
    126     points(i)=i;
    127125  }
    128126
    129   rl.fit(100, points);
     127  rl.fit(10, 100);
    130128
    131129  gslapi::vector y = rl.y_predicted();
Note: See TracChangeset for help on using the changeset viewer.