source: trunk/test/regression_test.cc @ 430

// $Id: regression_test.cc 430 2005-12-08 22:53:08Z peter $


#include <c++_tools/gslapi/matrix.h>
#include <c++_tools/statistics/KernelBox.h>
#include <c++_tools/statistics/Linear.h>
#include <c++_tools/statistics/LinearWeighted.h>
#include <c++_tools/statistics/Local.h>
#include <c++_tools/statistics/Naive.h>
#include <c++_tools/statistics/NaiveWeighted.h>
#include <c++_tools/statistics/Polynomial.h>

#include <cmath>

#include <fstream>
#include <iostream>


using namespace theplu;

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


int main(const int argc,const char* argv[])
{
  using namespace theplu;
  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)
  gslapi::vector x(4); x(0)=1970; x(1)=1980; x(2)=1990; x(3)=2000;
  gslapi::vector y(4); y(0)=12;   y(1)=11;   y(2)=14;   y(3)=13;
  gslapi::vector w(4); w(0)=0.1;  w(1)=0.2;  w(2)=0.3;  w(3)=0.4;

  *error << "testing regression::LinearWeighted" << std::endl;
  statistics::regression::LinearWeighted linear_w;
  linear_w.fit(x,y,w);
  double y_predicted=0;
  double y_predicted_err=0;
  linear_w.predict(1990,y_predicted,y_predicted_err);
  if (y_predicted!=12.8){
    *error << "regression_LinearWeighted: cannot reproduce fit." << std::endl;
    ok=false;
  }

  // testing regression::NaiveWeighted
  *error << "testing regression::NaiveWeighted" << std::endl;
  statistics::regression::NaiveWeighted naive_w;
  naive_w.fit(x,y,w);
  y_predicted=0;
  y_predicted_err=0;
  naive_w.predict(0.0,y_predicted,y_predicted_err);
  if (y_predicted!=(0.1*12+0.2*11+0.3*14+0.4*13)) {
    *error << "regression_NaiveWeighted: cannot reproduce fit." << std::endl;
    ok=false;
  }

  // testing regression::Local
  *error << "testing regression::Local" << std::endl;
  statistics::regression::KernelBox kb;
  statistics::regression::LinearWeighted rl;
  if (!Local_test(rl,kb)) {
    *error << "regression_Local: Linear cannot reproduce fit." << std::endl;
    ok=false;
  }
  statistics::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");
    gslapi::matrix data(s);
    gslapi::vector x(data.rows());
    gslapi::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));
    }

    statistics::regression::Polynomial polynomialfit(2);
    polynomialfit.fit(x,ln_y);
    gslapi::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;
  statistics::regression::Linear lin;
  
  *error << "testing regression::Naive" << std::endl;
  statistics::regression::Naive naive;

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

  return (ok ? 0 : -1);
}



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

  rl.fit(100, points);

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