Changeset 730 for trunk/yat/regression/LinearWeighted.cc

Timestamp:

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

Author:

Peter

Message:

fixes #167 and #160

File:

• trunk/yat/regression/LinearWeighted.cc (modified) (8 diffs)

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