Changeset 730

Timestamp:

Jan 6, 2007, 12:02:21 PM (16 years ago)

Author:

Peter

Message:

fixes #167 and #160

Location:

trunk

Files:

• test/regression_test.cc (modified) (6 diffs)
• yat/regression/Linear.cc (modified) (1 diff)
• yat/regression/LinearWeighted.cc (modified) (8 diffs)
• yat/regression/LinearWeighted.h (modified) (4 diffs)
• yat/regression/NaiveWeighted.cc (modified) (1 diff)

trunk/test/regression_test.cc

@@ -120 +120 @@
   *error << "  testing regression::LinearWeighted" << std::endl;
   regression::LinearWeighted linear_w;
+  ok = equal(linear, linear_w, error) && ok;
   linear_w.fit(x,y,w);
   double y_predicted = linear_w.predict(1990);
@@ -134 +135 @@
   regression::NaiveWeighted naive_w;
   regression::Naive naive;
-  ok = ok && equal(naive, naive_w, error);
+  ok = equal(naive, naive_w, error) && ok;
   naive_w.fit(x,y,w);
 
@@ -145 +146 @@
     ok=false;
   }
-
 
   // testing regression::Local
@@ -183 +183 @@
   }
 
-  *error << "  testing regression::Linear" << std::endl;
-  regression::Linear lin;
-  
   *error << "  testing regression::PolynomialWeighted" << std::endl;
   regression::PolynomialWeighted pol_weighted(2);
@@ -225 +222 @@
     *error << "Error: predict not equal" << std::endl;
   }
-  if (r.prediction_error2(2000) != rw.prediction_error2(2000)){
-    ok = false;
-    *error << "Error: prediction_error2 not equal non-weighted version." 
-           << std::endl;
-  }
-  if (r.r2() != rw.r2()){
+  if (fabs(r.prediction_error2(2000)-rw.prediction_error2(2000))>10e-7){
+    ok = false;
+    *error << "Error: prediction_error2 not equal non-weighted version.\n" 
+           << "       weighted: " << rw.prediction_error2(2000) << "\n"
+           << "   non-weighted: " << r.prediction_error2(2000)
+           << std::endl;
+  }
+  if (fabs(r.s2()-rw.s2())>10E-7){
     ok = false;
     *error << "Error: r2 not equal non-weighted version." << std::endl;
-  }
-  if (r.s2() != rw.s2()){
+    *error << "weighted r2 = " << rw.r2() << std::endl;
+    *error << "non-weighted r2 = " << r.r2() << std::endl;
+  }
+  if (fabs(r.s2()-rw.s2())>10E-7){
     ok = false;
     *error << "Error: s2 not equal non-weighted version." << std::endl;
-  }
-  if (r.standard_error2(2000) != rw.standard_error2(2000)){
+    *error << "weighted s2 = " << rw.s2() << std::endl;
+    *error << "non-weighted s2 = " << r.s2() << std::endl;
+  }
+  if (fabs(r.standard_error2(2000)-rw.standard_error2(2000))>10e-7){
     ok = false;
     *error << "Error: standard_error not equal non-weighted version." 
@@ -279 +282 @@
     *error << "Error: s2 not equal after rescaling.\n";
     *error << "       s2 = " << s2 << " and after doubling weights.\n";
-    *error << "       s2 = " << wr.s2() << "\n";
+    *error << "       s2 = " << wr.s2(2) << "\n";
   }
   if (wr.standard_error2(2000) != standard_error2){

trunk/yat/regression/Linear.cc

@@ -66 +66 @@
 
     alpha_ = ap_.y_averager().mean();
-    beta_ = ap_.covariance() / ap_.x_averager().variance();
+    beta_ = ap_.sum_xy_centered() / ap_.x_averager().sum_xx_centered();
 
     // calculating deviation between data and model

trunk/yat/regression/LinearWeighted.cc

@@ -26 +26 @@
 #include "yat/utility/vector.h"
 
+#include <cassert>
+
 namespace theplu {
 namespace yat {
@@ -32 +34 @@
   LinearWeighted::LinearWeighted(void)
     : OneDimensionalWeighted(), alpha_(0), alpha_var_(0), beta_(0),
-      beta_var_(0), m_x_(0), s2_(0)
+      beta_var_(0)
   {
   }
@@ -47 +49 @@
   double LinearWeighted::alpha_var(void) const
   {
-    return sqrt(alpha_var_);
+    return alpha_var_;
   }
 
@@ -57 +59 @@
   double LinearWeighted::beta_var(void) const
   {
-    return sqrt(beta_var_);
+    return beta_var_;
   }
 
@@ -64 +66 @@
                            const utility::vector& w)
   {
+    assert(x.size()==y.size());
+    assert(x.size()==w.size());
+
     // AveragerPairWeighted requires 2 weights but works only on the
     // product wx*wy, so we can send in w and a dummie to get what we
@@ -70 +75 @@
     ap_.add_values(x,y,utility::vector(x.size(),1),w);
 
-    // estimating the noise level. see attached document for motivation
-    // of the expression.
-    s2_= (syy()-sxy()*sxy()/sxx())/(w.sum()-2*(w*w)/w.sum()) ;
-    
     alpha_ = m_y();
     beta_ = sxy()/sxx();
-    alpha_var_ = ap_.y_averager().standard_error() * 
-      ap_.y_averager().standard_error();
-    beta_var_ = s2_/sxx();  
-    m_x_=m_x();
+
+    chisq_=0;
+    for (size_t i=0; i<x.size(); ++i){
+      double res = predict(x(i))-y(i);
+      chisq_ += w(i)*res*res;
+    }
+
+    alpha_var_ = s2()/ap_.y_averager().sum_w();
+    beta_var_ = s2()/sxx(); 
   }
 
@@ -92 +98 @@
   }
 
-  double LinearWeighted::prediction_error2(const double x, const double w) const
-  {
-    return sqrt(alpha_var_ + beta_var_*(x-m_x_)*(x-m_x_)+s2(w));
+  double LinearWeighted::predict(const double x) const
+  { 
+    return alpha_ + beta_ * (x-m_x()); 
   }
 
+  
   double LinearWeighted::s2(double w) const
   {
-    return s2_/w;
+    return chisq_/(w*(ap_.y_averager().n()-2));
   }
 
   double LinearWeighted::standard_error2(const double x) const
   {
-    return sqrt(alpha_var_ + beta_var_*(x-m_x_)*(x-m_x_) );
+    return alpha_var_ + beta_var_*(x-m_x())*(x-m_x());
   }
+
 
   double LinearWeighted::sxx(void) const
@@ -112 +120 @@
   }
 
+
   double LinearWeighted::sxy(void) const
   {
     return ap_.sum_xy_centered();
   }
+
 
   double LinearWeighted::syy(void) const

trunk/yat/regression/LinearWeighted.h

@@ -27 +27 @@
 #include "OneDimensionalWeighted.h"
 
-#include <cmath>
-
 namespace theplu {
 namespace yat {
@@ -55 +53 @@
     virtual ~LinearWeighted(void);
          
-    ///
-    /// @return the parameter \f$ \alpha \f$
-    ///
+    /**
+       \f$ alpha \f$ is estimated as \f$ \frac{\sum w_iy_i}{\sum w_i} \f$
+    
+       @return the parameter \f$ \alpha \f$
+    */
     double alpha(void) const;
 
-    ///
-    /// @return standard deviation of parameter \f$ \alpha \f$
-    ///
+    /**
+       Variance is estimated as \f$ \frac{s^2}{\sum w_i} \f$
+
+       @see s2()
+
+       @return variance of parameter \f$ \alpha \f$
+    */
     double alpha_var(void) const;
 
-    ///
-    /// @return the parameter \f$ \beta \f$
-    ///
+    /**
+       \f$ beta \f$ is estimated as \f$ \frac{\sum
+       w_i(y_i-m_y)(x_i-m_x)}{\sum w_i(x_i-m_x)^2} \f$
+    
+       @return the parameter \f$ \beta \f$
+    */
     double beta(void) const;
 
-    ///
-    /// @return standard deviation of parameter \f$ \beta \f$
-    ///
+    /**
+       Variance is estimated as \f$ \frac{s^2}{\sum w_i(x_i-m_x)^2} \f$
+
+       @see s2()
+
+       @return variance of parameter \f$ \beta \f$
+    */
     double beta_var(void) const;
     
@@ -91 +102 @@
     /// \f$ y =\alpha + \beta (x - m) \f$
     ///
-    double predict(const double x) const { return alpha_ + beta_ * (x-m_x_); }
-
-    ///
-    /// estimated squared deviation from predicted value for a new
-    /// data point in @a x with weight @a w
-    ///
-    double prediction_error2(const double x, const double w=1) const;
+    double predict(const double x) const;
 
     /**
@@ -126 +131 @@
     double beta_;
     double beta_var_;
-    double m_x_; // average of x values
-    double r2_; // coefficient of determination
-    double s2_;
-    double mse_;
   };
 

trunk/yat/regression/NaiveWeighted.cc

@@ -68 +68 @@
   double NaiveWeighted::standard_error2(const double x) const
   { 
-    return chisq_/((ap_.y_averager().n()-1)*ap_.sum_w());
+    return chisq_/( (ap_.y_averager().n()-1)*ap_.y_averager().sum_w() );
   }