Ignore:
Timestamp:
Dec 8, 2005, 5:20:36 PM (16 years ago)
Author:
Peter
Message:

improved doc in Score

File:
1 edited

Legend:

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

    r385 r428  
    3434  }
    3535
    36   void Linear::fit(const gslapi::vector& x,
    37                    const gslapi::vector& y,
    38                    const gslapi::vector& w)
    39   {
    40     double m_x = w*x /w.sum();
    41     double m_y = w*y /w.sum();
    42    
    43     double sxy = 0;
    44     for (size_t i=0; i<x.size(); i++)
    45       sxy += w(i)*(x(i)-m_x)*(y(i)-m_y);
    46 
    47     double sxx = 0;
    48     for (size_t i=0; i<x.size(); i++)
    49       sxx += w(i)*(x(i)-m_x)*(x(i)-m_x);
    50 
    51     double syy = 0;
    52     for (size_t i=0; i<y.size(); i++)
    53       syy += w(i)*(y(i)-m_y)*(y(i)-m_y);
    54 
    55     // estimating the noise level. see attached document for motivation
    56     // of the expression.
    57     s2_= (syy-sxy*sxy/sxx)/(w.sum()-2*(w*w)/w.sum()) ;
    58    
    59     alpha_ = m_y;
    60     beta_ = sxy/sxx;
    61     alpha_var_ = s2_/w.sum();
    62     beta_var_ = s2_/sxx; 
    63     m_x_=m_x;
     36  void Linear::predict(const double x, double& y, double& y_err)
     37  {
     38    y = alpha_ + beta_ * (x-m_x_);
     39    y_err = sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)+s2_/w );
     40    x_=x;
     41    y_=y;
     42    y_err_=y_err;
    6443  }
    6544
Note: See TracChangeset for help on using the changeset viewer.