Changeset 730 for trunk/yat

Timestamp:

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

Author:

Peter

Message:

fixes #167 and #160

Location:

trunk/yat/regression

Files:

• Linear.cc (modified) (1 diff)
• LinearWeighted.cc (modified) (8 diffs)
• LinearWeighted.h (modified) (4 diffs)
• NaiveWeighted.cc (modified) (1 diff)

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() );
   }