Ignore:
Timestamp:
Jan 6, 2007, 12:02:21 PM (15 years ago)
Author:
Peter
Message:

fixes #167 and #160

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/regression/LinearWeighted.cc

    r729 r730  
    2626#include "yat/utility/vector.h"
    2727
     28#include <cassert>
     29
    2830namespace theplu {
    2931namespace yat {
     
    3234  LinearWeighted::LinearWeighted(void)
    3335    : OneDimensionalWeighted(), alpha_(0), alpha_var_(0), beta_(0),
    34       beta_var_(0), m_x_(0), s2_(0)
     36      beta_var_(0)
    3537  {
    3638  }
     
    4749  double LinearWeighted::alpha_var(void) const
    4850  {
    49     return sqrt(alpha_var_);
     51    return alpha_var_;
    5052  }
    5153
     
    5759  double LinearWeighted::beta_var(void) const
    5860  {
    59     return sqrt(beta_var_);
     61    return beta_var_;
    6062  }
    6163
     
    6466                           const utility::vector& w)
    6567  {
     68    assert(x.size()==y.size());
     69    assert(x.size()==w.size());
     70
    6671    // AveragerPairWeighted requires 2 weights but works only on the
    6772    // product wx*wy, so we can send in w and a dummie to get what we
     
    7075    ap_.add_values(x,y,utility::vector(x.size(),1),w);
    7176
    72     // estimating the noise level. see attached document for motivation
    73     // of the expression.
    74     s2_= (syy()-sxy()*sxy()/sxx())/(w.sum()-2*(w*w)/w.sum()) ;
    75    
    7677    alpha_ = m_y();
    7778    beta_ = sxy()/sxx();
    78     alpha_var_ = ap_.y_averager().standard_error() *
    79       ap_.y_averager().standard_error();
    80     beta_var_ = s2_/sxx(); 
    81     m_x_=m_x();
     79
     80    chisq_=0;
     81    for (size_t i=0; i<x.size(); ++i){
     82      double res = predict(x(i))-y(i);
     83      chisq_ += w(i)*res*res;
     84    }
     85
     86    alpha_var_ = s2()/ap_.y_averager().sum_w();
     87    beta_var_ = s2()/sxx();
    8288  }
    8389
     
    9298  }
    9399
    94   double LinearWeighted::prediction_error2(const double x, const double w) const
    95   {
    96     return sqrt(alpha_var_ + beta_var_*(x-m_x_)*(x-m_x_)+s2(w));
     100  double LinearWeighted::predict(const double x) const
     101  { 
     102    return alpha_ + beta_ * (x-m_x());
    97103  }
    98104
     105 
    99106  double LinearWeighted::s2(double w) const
    100107  {
    101     return s2_/w;
     108    return chisq_/(w*(ap_.y_averager().n()-2));
    102109  }
    103110
    104111  double LinearWeighted::standard_error2(const double x) const
    105112  {
    106     return sqrt(alpha_var_ + beta_var_*(x-m_x_)*(x-m_x_) );
     113    return alpha_var_ + beta_var_*(x-m_x())*(x-m_x());
    107114  }
     115
    108116
    109117  double LinearWeighted::sxx(void) const
     
    112120  }
    113121
     122
    114123  double LinearWeighted::sxy(void) const
    115124  {
    116125    return ap_.sum_xy_centered();
    117126  }
     127
    118128
    119129  double LinearWeighted::syy(void) const
Note: See TracChangeset for help on using the changeset viewer.