source: trunk/test/regression_test.cc @ 741

// $Id: regression_test.cc 741 2007-01-13 14:41:40Z peter $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "yat/regression/KernelBox.h"
#include "yat/regression/Linear.h"
#include "yat/regression/LinearWeighted.h"
#include "yat/regression/Local.h"
#include "yat/regression/Naive.h"
#include "yat/regression/NaiveWeighted.h"
#include "yat/regression/Polynomial.h"
#include "yat/regression/PolynomialWeighted.h"
#include "yat/utility/matrix.h"
#include "yat/utility/vector.h"

#include <cmath>

#include <fstream>
#include <iostream>


using namespace theplu::yat;

bool equal(regression::OneDimensional&, regression::OneDimensionalWeighted&, 
           std::ostream*); 

bool multidim(std::ostream* error);

bool unity_weights(regression::OneDimensional&, 
                   regression::OneDimensionalWeighted&, 
                   const utility::vector&, const utility::vector&,
                   std::ostream*); 

bool rescale_weights(regression::OneDimensionalWeighted&, 
                     const utility::vector&, const utility::vector&,
                     std::ostream*); 

bool zero_weights(regression::OneDimensionalWeighted&, 
                  const utility::vector&, const utility::vector&,
                  std::ostream*); 


bool Local_test(regression::OneDimensionalWeighted&, 
                regression::Kernel&);

int main(const int argc,const char* argv[])
{
  std::ostream* error;
  if (argc>1 && argv[1]==std::string("-v"))
    error = &std::cerr;
  else {
    error = new std::ofstream("/dev/null");
    if (argc>1)
      std::cout << "regression_test -v : for printing extra information\n";
  }
  *error << "  testing regression" << std::endl;
  bool ok = true;

  // test data for Linear and Naive (Weighted and non-weighted)
  utility::vector x(4); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000;
  utility::vector y(4); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;
  utility::vector w(4); w(0)=0.1;  w(1)=0.2;  w(2)=0.3;  w(3)=0.4;

  // Comparing linear and polynomial(1)
  regression::Linear linear;
  linear.fit(x,y);
  regression::Polynomial polynomial(1);
  polynomial.fit(x,y);
  if ( fabs(linear.beta()-polynomial.fit_parameters()(1))>0.0001 ){
    *error << "error: beta and fit_parameters(1) not equal" << std::endl;
    *error << "       beta = " << linear.beta() << std::endl;
    *error << "       fit_parameters(1) = " 
           << polynomial.fit_parameters()(1) << std::endl;
    ok = false;
  }
  if ( fabs(polynomial.fit_parameters()(0)-linear.alpha()+
            linear.beta()*1985)>0.0001){
    *error << "error: fit_parameters(0) = " 
           << polynomial.fit_parameters()(0)<< std::endl;
    *error << "error: alpha-beta*m_x = " 
           << linear.alpha()-linear.beta()*1985 << std::endl;
    ok = false;
  }
  if ( fabs(polynomial.chisq()-linear.chisq())>0.0001){
    *error << "error: chisq not same in linear and polynomial(1)" 
           << std::endl;
    ok = false;
  }
  if ( fabs(polynomial.predict(1.0)-linear.predict(1.0))>0.0001){
    *error << "error: predict not same in linear and polynomial(1)" 
           << std::endl;
    ok = false;
  }
  if ( fabs(polynomial.standard_error2(1985)-linear.standard_error2(1985))
       >0.0001){
    *error << "error: standard_error not same in linear and polynomial(1)" 
           << "\n  polynomial: " << polynomial.standard_error2(1985)
           << "\n  linear: " << linear.standard_error2(1985)
           << std::endl;
    ok = false;
  }

  *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);
  double y_predicted_err = linear_w.prediction_error2(1990);
  if (fabs(y_predicted-12.8)>0.001){
    *error << "error: cannot reproduce fit." << std::endl;
    *error << "predicted value: " << y_predicted << " expected 12.8" 
           << std::endl;
    ok=false;
  }

  // testing regression::NaiveWeighted
  *error << "  testing regression::NaiveWeighted" << std::endl;
  regression::NaiveWeighted naive_w;
  regression::Naive naive;
  ok = equal(naive, naive_w, error) && ok;
  naive_w.fit(x,y,w);

  y_predicted=naive_w.predict(0.0);
  y_predicted_err=naive_w.prediction_error2(0.0);
  if (y_predicted!=(0.1*12+0.2*11+0.3*14+0.4*13)) {
    *error << "regression_NaiveWeighted: cannot reproduce fit." << std::endl;
    *error << "returned value: " << y_predicted << std::endl;
    *error << "expected: " << 0.1*12+0.2*11+0.3*14+0.4*13 << std::endl;
    ok=false;
  }

  // testing regression::Local
  *error << "  testing regression::Local" << std::endl;
  regression::KernelBox kb;
  regression::LinearWeighted rl;
  if (!Local_test(rl,kb)) {
    *error << "regression_Local: Linear cannot reproduce fit." << std::endl;
    ok=false;
  }
  regression::NaiveWeighted rn;
  if (!Local_test(rn,kb)) {
    *error << "regression_Local: Naive cannot reproduce fit." << std::endl;
    ok=false;
  }

  // testing regression::Polynomial
  *error << "  testing regression::Polynomial" << std::endl;
  {
    std::ifstream s("data/regression_gauss.data");
    utility::matrix data(s);
    utility::vector x(data.rows());
    utility::vector ln_y(data.rows());
    for (size_t i=0; i<data.rows(); ++i) {
      x(i)=data(i,0);
      ln_y(i)=log(data(i,1));
    }

    regression::Polynomial polynomialfit(2);
    polynomialfit.fit(x,ln_y);
    utility::vector fit=polynomialfit.fit_parameters();
    if (fabs(fit[0]-1.012229646706 + fit[1]-0.012561322528 +
             fit[2]+1.159674470130)>1e-11) {  // Jari, fix number!
      *error << "regression_Polynomial: cannot reproduce fit." << std::endl;
      ok=false;
    }
  }

  *error << "  testing regression::PolynomialWeighted" << std::endl;
  regression::PolynomialWeighted pol_weighted(2);
  regression::Polynomial polynomial2(2);
  ok = equal(polynomial2, pol_weighted, error) && ok;

  ok = multidim(error) && ok;

  if (!ok)
    *error << "Test failed\n" << std::endl;

  if (error!=&std::cerr)
    delete error;

  return (ok ? 0 : -1);
}


bool equal(regression::OneDimensional& r, 
           regression::OneDimensionalWeighted& wr, 
           std::ostream* error)
{
  bool ok=true;
  utility::vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
  utility::vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;

  ok = unity_weights(r, wr, x, y, error) && ok;
  ok = rescale_weights(wr, x, y, error) && ok;
  ok = zero_weights(wr, x, y, error) && ok;
  return ok;
}


bool unity_weights(regression::OneDimensional& r, 
                   regression::OneDimensionalWeighted& rw, 
                   const utility::vector& x, const utility::vector& y,
                   std::ostream* error)
{
  *error << "    testing unity weights equal to non-weighted version.\n"; 
  bool ok=true;
  utility::vector w(x.size(), 1.0);
  r.fit(x,y);
  rw.fit(x,y,w);
  if (r.predict(2000) != rw.predict(2000)){
    ok = false;
    *error << "Error: predict not equal\n" 
           << "   weighted: " << rw.predict(2000) << "\n"
           << "   non-weighted: " << r.predict(2000)
           << std::endl;
  }
  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;
    *error << "weighted r2 = " << rw.r2() << std::endl;
    *error << "non-weighted r2 = " << r.r2() << std::endl;
  }
  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;
  }
  return ok;
} 


bool rescale_weights(regression::OneDimensionalWeighted& wr, 
                     const utility::vector& x, const utility::vector& y,
                     std::ostream* error)
{
  *error << "    testing rescaling weights.\n"; 
  bool ok = true;
  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
  wr.fit(x,y,w);
  double predict = wr.predict(2000);
  double prediction_error2 = wr.prediction_error2(2000);
  double r2 = wr.r2();
  double s2 = wr.s2();
  double standard_error2 = wr.standard_error2(2000);

  w.scale(2);
  wr.fit(x,y,w);
  if (fabs(wr.predict(2000)-predict)>10e-11){
    ok = false;
    *error << "Error: predict not equal after rescaling.\n";
    *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";
    *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; 
}


bool zero_weights(regression::OneDimensionalWeighted& wr, 
                  const utility::vector& x, const utility::vector& y,
                  std::ostream* error)
{
  *error << "    testing zero weights equal to missing value.\n"; 
  bool ok = true;
  utility::vector w(5);  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0;
  wr.fit(x,y,w);
  double predict = wr.predict(2000);
  double prediction_error2 = wr.prediction_error2(2000);
  double r2 = wr.r2();
  double s2 = wr.s2();
  double standard_error2 = wr.standard_error2(2000);

  utility::vector x2(4);
  utility::vector y2(4);
  utility::vector w2(4);
  for (size_t i=0; i<4; ++i){
    x2(i) = x(i);
    y2(i) = y(i);
    w2(i) = w(i);
  }

  wr.fit(x2,y2,w2);
  if (wr.predict(2000) != predict){
    ok = false;
    *error << "Error: predict not equal.\n";
  }
  if (wr.prediction_error2(2000) != prediction_error2){
    ok = false;
    *error << "Error: prediction_error2 not equal.\n";
  }
  if (wr.r2() != r2){
    ok = false;
    *error << "Error: r2 not equal.\n";
    *error << "   r2: " << r2 << "\n";
    *error << "   r2: " << wr.r2() << "\n";
  }
  if (wr.s2() != s2){
    ok = false;
    *error << "Error: s2 not equal.\n";
  }
  if (wr.standard_error2(2000) != standard_error2){
    ok = false;
    *error << "Error: standard_error2 not equal.\n";
  }
  return ok; 
}


bool multidim(std::ostream* error)
{
  bool ok = true;
  *error << "  testing regression::MultiDimensionalWeighted" << std::endl;
  utility::vector x(5); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000; x(4)=2010;
  utility::vector y(5); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;   y(4)=15;
  utility::vector w(5,1.0);
  
  utility::matrix data(5,3);
  for (size_t i=0; i<data.rows(); ++i){
    data(i,0)=1;
    data(i,1)=x(i);
    data(i,2)=x(i)*x(i);
  }
  regression::MultiDimensional md;
  md.fit(data,y);
  regression::MultiDimensionalWeighted mdw;
  mdw.fit(data,y,w);
  utility::vector z(3,1);
  z(1)=2000;
  z(2)=2000*2000;
  if (md.predict(z) != mdw.predict(z)){
    ok = false;
    *error << "Error: predict not equal\n" 
           << "   weighted: " << mdw.predict(z) << "\n"
           << "   non-weighted: " << md.predict(z)
           << std::endl;
  }

  if (md.standard_error2(z) != mdw.standard_error2(z)){
    ok = false;
    *error << "Error: standard_error2 not equal\n" 
           << "   weighted: " << mdw.standard_error2(z) << "\n"
           << "   non-weighted: " << md.standard_error2(z)
           << std::endl;
  }
  if (fabs(md.prediction_error2(z)-mdw.prediction_error2(z,1.0))>10e-7){
    ok = false;
    *error << "Error: prediction_error2 not equal\n" 
           << "   weighted: " << mdw.prediction_error2(z,1.0) << "\n"
           << "   non-weighted: " << md.prediction_error2(z)
           << std::endl;
  }

  w(0)=1.0;  w(1)=1.0;  w(2)=0.5;  w(3)=0.2;  w(4)=0.2;
  mdw.fit(data,y,w);
  double predict = mdw.predict(z);
  double prediction_error2 = mdw.prediction_error2(z, 1.0);
  double s2 = mdw.s2(1.0);
  double standard_error2 = mdw.standard_error2(z);

  w.scale(2);
  mdw.fit(data,y,w);
  if (fabs(mdw.predict(z)-predict)>10e-10){
    ok = false;
    *error << "Error: predict not equal after rescaling.\n";
    *error << "   predict = " << predict << " and after doubling weights.\n";
    *error << "   predict = " << mdw.predict(z) << "\n";
  }
  if (fabs(mdw.prediction_error2(z,2)-prediction_error2)>10e-7){ 
    ok = false;
    *error << "Error: prediction_error2 not equal after rescaling.\n";
    *error << "       predict_error2 = " << prediction_error2 
           << " and after doubling weights.\n";
    *error << "       predict_error2 = " << mdw.prediction_error2(z,2) << "\n";
  }
  if (fabs(mdw.s2(2)-s2)>10e-10){
    ok = false;
    *error << "Error: s2 not equal after rescaling.\n";
    *error << "       s2 = " << s2 << " and after doubling weights.\n";
    *error << "       s2 = " << mdw.s2(2) << "\n";
  }
  if (fabs(mdw.standard_error2(z)-standard_error2)>10e-7){
    ok = false;
    *error << "Error: standard_error2 not equal after rescaling.\n";
    *error << " standard_error2 = " << standard_error2 
           << " and after doubling weights.\n";
    *error << " standard_error2 = " << mdw.standard_error2(z) << "\n";
  }

  return ok;
}


bool Local_test(regression::OneDimensionalWeighted& r, 
                regression::Kernel& k)
{
  regression::Local rl(r,k);
  for (size_t i=0; i<1000; i++){
    rl.add(i, 10);
  }

  rl.fit(10, 100);

  utility::vector y = rl.y_predicted();
  for (size_t i=0; i<y.size(); i++) 
    if (y(i)!=10.0){
      return false;
    }
  return true;
}