Changeset 730 for trunk/yat/regression


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

fixes #167 and #160

Location:
trunk/yat/regression
Files:
4 edited

Legend:

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

    r729 r730  
    6666
    6767    alpha_ = ap_.y_averager().mean();
    68     beta_ = ap_.covariance() / ap_.x_averager().variance();
     68    beta_ = ap_.sum_xy_centered() / ap_.x_averager().sum_xx_centered();
    6969
    7070    // calculating deviation between data and model
  • 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
  • trunk/yat/regression/LinearWeighted.h

    r729 r730  
    2727#include "OneDimensionalWeighted.h"
    2828
    29 #include <cmath>
    30 
    3129namespace theplu {
    3230namespace yat {
     
    5553    virtual ~LinearWeighted(void);
    5654         
    57     ///
    58     /// @return the parameter \f$ \alpha \f$
    59     ///
     55    /**
     56       \f$ alpha \f$ is estimated as \f$ \frac{\sum w_iy_i}{\sum w_i} \f$
     57   
     58       @return the parameter \f$ \alpha \f$
     59    */
    6060    double alpha(void) const;
    6161
    62     ///
    63     /// @return standard deviation of parameter \f$ \alpha \f$
    64     ///
     62    /**
     63       Variance is estimated as \f$ \frac{s^2}{\sum w_i} \f$
     64
     65       @see s2()
     66
     67       @return variance of parameter \f$ \alpha \f$
     68    */
    6569    double alpha_var(void) const;
    6670
    67     ///
    68     /// @return the parameter \f$ \beta \f$
    69     ///
     71    /**
     72       \f$ beta \f$ is estimated as \f$ \frac{\sum
     73       w_i(y_i-m_y)(x_i-m_x)}{\sum w_i(x_i-m_x)^2} \f$
     74   
     75       @return the parameter \f$ \beta \f$
     76    */
    7077    double beta(void) const;
    7178
    72     ///
    73     /// @return standard deviation of parameter \f$ \beta \f$
    74     ///
     79    /**
     80       Variance is estimated as \f$ \frac{s^2}{\sum w_i(x_i-m_x)^2} \f$
     81
     82       @see s2()
     83
     84       @return variance of parameter \f$ \beta \f$
     85    */
    7586    double beta_var(void) const;
    7687   
     
    91102    /// \f$ y =\alpha + \beta (x - m) \f$
    92103    ///
    93     double predict(const double x) const { return alpha_ + beta_ * (x-m_x_); }
    94 
    95     ///
    96     /// estimated squared deviation from predicted value for a new
    97     /// data point in @a x with weight @a w
    98     ///
    99     double prediction_error2(const double x, const double w=1) const;
     104    double predict(const double x) const;
    100105
    101106    /**
     
    126131    double beta_;
    127132    double beta_var_;
    128     double m_x_; // average of x values
    129     double r2_; // coefficient of determination
    130     double s2_;
    131     double mse_;
    132133  };
    133134
  • trunk/yat/regression/NaiveWeighted.cc

    r729 r730  
    6868  double NaiveWeighted::standard_error2(const double x) const
    6969  {
    70     return chisq_/((ap_.y_averager().n()-1)*ap_.sum_w());
     70    return chisq_/( (ap_.y_averager().n()-1)*ap_.y_averager().sum_w() );
    7171  }
    7272
Note: See TracChangeset for help on using the changeset viewer.