source: trunk/test/regression_test.cc @ 726

// $Id: regression_test.cc 726 2007-01-04 14:38:56Z 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 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_error(1985)-linear.standard_error(1985))
       >0.0001){
    *error << "error: standard_error not same in linear and polynomial(1)" 
           << "\n  polynomial: " << polynomial.standard_error(1.0)
           << "\n  linear: " << linear.standard_error(1.0)
           << "\n  alpha_var: " << linear.alpha_var()
           << "\n  beta_var: " << linear.beta_var()
           << "\n  covariance: " << polynomial.covariance()(0,0)
           << " " << polynomial.covariance()(0,1) << "\n"
           << " " << polynomial.covariance()(1,0) 
           << " " << polynomial.covariance()(1,1) 
           << std::endl;
    ok = false;
  }

  *error << "testing regression::LinearWeighted" << std::endl;
  regression::LinearWeighted linear_w;
  linear_w.fit(x,y,w);
  double y_predicted = linear_w.predict(1990);
  double y_predicted_err = linear_w.prediction_error(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;
  naive_w.fit(x,y,w);

  y_predicted=naive_w.predict(0.0);
  y_predicted_err=naive_w.prediction_error(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::Linear" << std::endl;
  regression::Linear lin;
  
  *error << "testing regression::Naive" << std::endl;
  regression::Naive naive;

  *error << "testing regression::Polynomial" << std::endl;
  regression::Polynomial pol(2);

  *error << "testing regression::PolynomialWeighted" << std::endl;
  regression::PolynomialWeighted pol_weighted(2);

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

  return (ok ? 0 : -1);
}



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