Changeset 741

Timestamp:

Jan 13, 2007, 3:41:40 PM (17 years ago)

Author:

Peter

Message:

fixes #161 and #164

Location:

trunk

Files:

• test/regression_test.cc (modified) (6 diffs)
• yat/regression/MultiDimensionalWeighted.cc (modified) (1 diff)
• yat/regression/MultiDimensionalWeighted.h (modified) (1 diff)
• yat/regression/Polynomial.cc (modified) (1 diff)
• yat/regression/PolynomialWeighted.cc (modified) (3 diffs)
• yat/regression/PolynomialWeighted.h (modified) (1 diff)

trunk/test/regression_test.cc

@@ -188 +188 @@
   regression::PolynomialWeighted pol_weighted(2);
   regression::Polynomial polynomial2(2);
-  // ok = equal(polynomial2, pol_weighted, error) && ok;
+  ok = equal(polynomial2, pol_weighted, error) && ok;
 
   ok = multidim(error) && ok;
@@ -234 +234 @@
            << std::endl;
   }
-  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){
+  if (fabs(r.s2()-rw.s2(1.0))>10E-7){
+    ok = false;
+    *error << "Error: s2 not equal non-weighted version." << std::endl;
+    *error << "weighted s2 = " << rw.s2(1.0) << 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." 
+           << std::endl;
+  }
+  if (fabs(r.r2()-rw.r2())>10E-7){
     ok = false;
     *error << "Error: r2 not equal non-weighted version." << std::endl;
@@ -247 +251 @@
     *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;
-    *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." 
+  if (fabs(r.prediction_error2(2000)-rw.prediction_error2(2000,1))>10e-7){
+    ok = false;
+    *error << "Error: prediction_error2 not equal non-weighted version.\n" 
+           << "       weighted: " << rw.prediction_error2(2000,1) << "\n"
+           << "   non-weighted: " << r.prediction_error2(2000)
            << std::endl;
   }
@@ -278 +278 @@
   w.scale(2);
   wr.fit(x,y,w);
-  if (wr.predict(2000) != predict){
+  if (fabs(wr.predict(2000)-predict)>10e-11){
     ok = false;
     *error << "Error: predict not equal after rescaling.\n";
-  }
-  if (wr.prediction_error2(2000,2) != prediction_error2){
-    ok = false;
-    *error << "Error: prediction_error2 not equal after rescaling.\n";
-  }
-  if (wr.r2() != r2){
-    ok = false;
-    *error << "Error: r2 not equal after rescaling.\n";
-  }
-  if (wr.s2(2) != s2){
+    *error << "       predict = " << predict 
+           << " and after doubling weights.\n";
+    *error << "       predict = " << wr.predict(2000) << "\n";
+  }
+  if (fabs(wr.s2(2)-s2)>10e-11){
     ok = false;
     *error << "Error: s2 not equal after rescaling.\n";
     *error << "       s2 = " << s2 << " and after doubling weights.\n";
     *error << "       s2 = " << wr.s2(2) << "\n";
-  }
-  if (wr.standard_error2(2000) != standard_error2){
+    *error << "difference " << s2-wr.s2(2.0) << std::endl;
+  }
+  if (fabs(wr.standard_error2(2000)-standard_error2)>10e-6){
     ok = false;
     *error << "Error: standard_error2 not equal after rescaling.\n";
+    *error << "       standard_error2 = " << standard_error2 
+           << " and after doubling weights.\n";
+    *error << "       standard_error2 = " 
+           << wr.standard_error2(2000) << "\n";
+    *error << " difference: " << wr.standard_error2(2000)-standard_error2
+           << std::endl;
+  }
+  if (fabs(wr.r2()-r2)>10e-6){
+    ok = false;
+    *error << "Error: r2 not equal after rescaling.\n";
+  }
+  if (fabs(wr.prediction_error2(2000,2)-prediction_error2)>10e-6){
+    ok = false;
+    *error << "Error: prediction_error2 not equal after rescaling.\n";
+    *error << "       prediction_error2 = " << prediction_error2 
+           << " and after doubling weights.\n";
+    *error << "       prediction_error2 = " 
+           << wr.prediction_error2(2000,2) << "\n";
   }
   return ok; 
@@ -339 +353 @@
     ok = false;
     *error << "Error: r2 not equal.\n";
+    *error << "   r2: " << r2 << "\n";
+    *error << "   r2: " << wr.r2() << "\n";
   }
   if (wr.s2() != s2){
@@ -411 +427 @@
     *error << "   predict = " << mdw.predict(z) << "\n";
   }
-  if (fabs(mdw.prediction_error2(z,2)-prediction_error2)>10e-7){
+  if (fabs(mdw.prediction_error2(z,2)-prediction_error2)>10e-7){ 
     ok = false;
     *error << "Error: prediction_error2 not equal after rescaling.\n";

trunk/yat/regression/MultiDimensionalWeighted.cc

@@ -43 +43 @@
       gsl_multifit_linear_free(work_);
   }
+
+
+  double MultiDimensionalWeighted::chisq() const
+  {
+    return chisquare_;
+  }
+
 
   void MultiDimensionalWeighted::fit(const utility::matrix& x, 

trunk/yat/regression/MultiDimensionalWeighted.h

@@ -52 +52 @@
 
     ///
+    /// @return sum of squared residuals
+    ///
+    double chisq(void) const;
+
+    ///
     /// @see gsl_multifit_wlinear
     ///

trunk/yat/regression/Polynomial.cc

@@ -49 +49 @@
   void Polynomial::fit(const utility::vector& x, const utility::vector& y)
   {
+    ap_.add_values(x,y);
     utility::matrix X=utility::matrix(x.size(),power_+1,1);
     for (size_t i=0; i<X.rows(); ++i)

trunk/yat/regression/PolynomialWeighted.cc

@@ -47 +47 @@
     assert(x.size()==y.size());
     assert(y.size()==w.size());
+    ap_.reset();
+    // 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
+    // want.
+    ap_.add_values(x,y,utility::vector(x.size(),1),w);
     utility::matrix X=utility::matrix(x.size(),power_+1,1);
     for (size_t i=0; i<X.rows(); ++i)
@@ -52 +57 @@
         X(i,j)=X(i,j-1)*x(i);
     md_.fit(X,y,w);
+    chisq_=md_.chisq();
   }
 
@@ -74 +80 @@
   }
 
-  double PolynomialWeighted::prediction_error2(const double x, 
-                                               const double w) const
-  {
-    utility::vector vec(power_+1,1);
-    for (size_t i=1; i<=power_; ++i)
-      vec(i) = vec(i-1)*x;
-    return md_.prediction_error2(vec, w);
-  }
-
   double PolynomialWeighted::standard_error2(const double x) const
   {

trunk/yat/regression/PolynomialWeighted.h

@@ -85 +85 @@
 
     ///
-    /// @return expected prediction error for a new data point in @a x
-    /// with weight @a w
-    ///
-    double prediction_error2(const double x, const double w=1) const;
-
-    ///
     /// @return error of model value in @a x
     ///