Changeset 726 for trunk/yat

Timestamp:

Jan 4, 2007, 3:38:56 PM (16 years ago)

Author:

Peter

Message:

fixes #165 added test checking Linear Regression is equivalent to Polynomial regression of degree one.

Location:

trunk/yat/regression

Files:

• Linear.cc (modified) (5 diffs)
• Linear.h (modified) (4 diffs)
• Local.h (modified) (1 diff)
• MultiDimensional.cc (modified) (1 diff)
• MultiDimensional.h (modified) (1 diff)
• Naive.cc (modified) (1 diff)
• Naive.h (modified) (1 diff)
• OneDimensional.h (modified) (1 diff)
• Polynomial.cc (modified) (1 diff)
• Polynomial.h (modified) (2 diffs)

trunk/yat/regression/Linear.cc

@@ -32 +32 @@
   Linear::Linear(void)
     : OneDimensional(), alpha_(0), alpha_var_(0), beta_(0), beta_var_(0),
-      chisq_(0), m_x_(0)
+      chisq_(0)
   {
   }
@@ -45 +45 @@
   }
 
-  double Linear::alpha_err(void) const
+  double Linear::alpha_var(void) const
   {
-    return sqrt(alpha_var_);
+    return alpha_var_;
   }
 
@@ -55 +55 @@
   }
 
-  double Linear::beta_err(void) const
+  double Linear::beta_var(void) const
   {
-    return sqrt(beta_var_);
+    return beta_var_;
   }
 
@@ -75 +75 @@
 
     // calculating deviation between data and model
-    chisq_ = ( (ap_.y_averager().sum_xx_centered() - ap_.sum_xy_centered()*
-              ap_.sum_xy_centered()/ap_.x_averager().sum_xx_centered() ) / 
-               (x.size()-2) );
-    r2_= 1-chisq_/ap_.x_averager().variance();
-    alpha_var_ = chisq_ / x.size();
-    beta_var_ = chisq_ / ap_.x_averager().sum_xx_centered();
-    m_x_ = ap_.x_averager().mean();
+    chisq_ = (ap_.y_averager().sum_xx_centered() - ap_.sum_xy_centered()*
+              ap_.sum_xy_centered()/ap_.x_averager().sum_xx_centered() );
+    r2_= 1-chisq_/ap_.x_averager().sum_xx_centered();
+    alpha_var_ = s2() / x.size();
+    beta_var_ = s2() / ap_.x_averager().sum_xx_centered();
   }
 
   double Linear::predict(const double x) const
   { 
-    return alpha_ + beta_ * (x-m_x_); 
+    return alpha_ + beta_ * (x - ap_.x_averager().mean()); 
   }
 
@@ -94 +92 @@
   }
 
+  double Linear::s2(void) const
+  {
+    return chisq()/(ap_.n()-2);
+  }
+
   double Linear::standard_error(const double x) const
   {
-    return sqrt( alpha_var_+beta_var_*(x-m_x_)*(x-m_x_)); 
+    return sqrt( alpha_var_+beta_var_*(x-ap_.x_averager().mean())*
+                 (x-ap_.x_averager().mean()) ); 
   }
 

trunk/yat/regression/Linear.h

@@ -65 +65 @@
 
     /**
-       The standard deviation is estimated as \f$
-       \sqrt{\frac{\chi^2}{n}} \f$
+       The standard deviation is estimated as \f$ \sqrt{\frac{s^2}{n}}
+       \f$ where \f$ s^2 = \frac{\sum \epsilon^2}{n-2} \f$
        
        @return standard deviation of parameter \f$ \alpha \f$
     */
-    double alpha_err(void) const;
+    double alpha_var(void) const;
 
     /**
@@ -81 +81 @@
 
     /**
-       The standard deviation is estimated as \f$
-       \sqrt{\frac{\chi^2}{\sum (x_i-m_x)^2}} \f$
-       
+       The standard deviation is estimated as \f$ \frac{s^2}{\sum
+       (x-m_x)^2} \f$ where \f$ s^2 = \frac{\sum \epsilon^2}{n-2} \f$
+
        @return standard deviation of parameter \f$ \beta \f$
     */
-    double beta_err(void) const;
+    double beta_var(void) const;
 
   /**
-       @brief Mean Squared Error
-    
-       Chisq is calculated as \f$ \frac{\sum
-       (y_i-\alpha-\beta(x_i-m_x))^2}{n-2} \f$
+       Chi-squared is calculated as \f$ \sum
+       (y_i-\alpha-\beta(x_i-m_x))^2 \f$
     */
     double chisq(void) const;
@@ -129 +127 @@
     Linear(const Linear&);
 
+    double s2(void) const;
+
     double alpha_;
     double alpha_var_;
@@ -134 +134 @@
     double beta_var_;
     double chisq_;
-    double m_x_; // average of x values
     double r2_; // coefficient of determination
   };

trunk/yat/regression/Local.h

@@ -102 +102 @@
 
   ///
-  /// The output operator for the RegressionLocal class.
+  /// The output operator for the Regression::Local class.
   ///
   std::ostream& operator<<(std::ostream&, const Local& );

trunk/yat/regression/MultiDimensional.cc

@@ -41 +41 @@
     if (work_)
       gsl_multifit_linear_free(work_);
+  }
+
+
+  const utility::matrix& MultiDimensional::covariance(void) const
+  {
+    return covariance_;
   }
 

trunk/yat/regression/MultiDimensional.h

@@ -52 +52 @@
 
     ///
+    /// @brief covariance of parameters
+    ///
+    const utility::matrix& covariance(void) const;
+
+    ///
     /// Function fitting parameters of the linear model by miminizing
     /// the quadratic deviation between model and data.

trunk/yat/regression/Naive.cc

@@ -47 +47 @@
   double Naive::chisq(void) const 
   { 
-    return ap_.y_averager().sum_xx_centered()/(ap_.n()-1); 
+    return ap_.y_averager().sum_xx_centered(); 
   }
   

trunk/yat/regression/Naive.h

@@ -58 +58 @@
          
     /** 
-        \f$\frac{1}{N-1} \sum (x_i-m)^2 \f$
-    
-        @brief Mean Squared Error
+        Chi-squared \f$ \sum (x_i-m)^2 \f$
     */
     double chisq(void) const; 

trunk/yat/regression/OneDimensional.h

@@ -56 +56 @@
  
     /**
-       @brief Mean Squared Error
+       @brief Chi-squared
        
-       Mean Squared Error is defined as the \f$ \frac
-       {\sum{(\hat{y_i}-y_i)^2}}{df} \f$ where \f$ df \f$ number of
-       degree of freedom typically is number data points minus number
-       of paramers in model.
+       Chi-squared is defined as the \f$ \frac
+       {\sum{(\hat{y_i}-y_i)^2}}{1} \f$
     */
     virtual double chisq(void) const=0;

trunk/yat/regression/Polynomial.cc

@@ -47 +47 @@
 
 
+  const utility::matrix& Polynomial::covariance(void) const
+  {
+    return md_.covariance();
+  }
+
+
   void Polynomial::fit(const utility::vector& x, const utility::vector& y)
   {

trunk/yat/regression/Polynomial.h

@@ -58 +58 @@
 
     ///
+    /// @brief covariance of parameters
+    ///
+    const utility::matrix& covariance(void) const;
+
+    ///
     /// Fit the model by minimizing the mean squared deviation between
     /// model and data.
@@ -71 +76 @@
 
     ///
-    /// @brief Variance of residuals
+    /// @brief Sum of squared residuals
     ///
     double chisq(void) const;