source: trunk/test/test_regression_linear.cc @ 212

// $Id: test_regression_linear.cc 212 2004-11-03 15:07:00Z jari $

// C++ tools include
////////////////////
#include "RegressionLinear.h"
#include "vector.h"


// Standard includes
////////////////////
#include <gsl/gsl_fit.h>
#include <iostream>


using namespace std;

int main()

{  
  bool ok=true;

  theplu::gslapi::vector x(4);
  x(0)=1970;
  x(1)=1980;
  x(2)=1990;
  x(3)=2000;
  theplu::gslapi::vector y(4);
  y(0)=12;
  y(1)=11;
  y(2)=14;
  y(3)=13;
  theplu::gslapi::vector w(4);
  w(0)=0.1;
  w(1)=0.2;
  w(2)=0.3;
  w(3)=0.4;

  double c0, c1, cov00, cov01, cov11, chisq;
  int n=4;
  int i=4;
  gsl_fit_wlinear(x.gsl_vector_pointer()->data, x.gsl_vector_pointer()->stride,
                  w.gsl_vector_pointer()->data, w.gsl_vector_pointer()->stride,
                  y.gsl_vector_pointer()->data, y.gsl_vector_pointer()->stride,
                  x.gsl_vector_pointer()->size,
                  &c0, &c1,
                  &cov00, &cov01, &cov11, &chisq);

  printf("# best fit: Y = %g + %g X\n", c0, c1);
  printf("# covariance matrix:\n");
  printf("# [ %g, %g\n#   %g, %g]\n", cov00, cov01, cov01, cov11);
  printf("# chisq = %g\n", chisq);

  for (i = 0; i < n; i++)
    printf("data: %g %g %g\n", x[i], y[i], 1/sqrt(w[i]));

  printf("\n");

  for (i = -30 ; i < 130 ; i++)
    {
      double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
      double yf, yf_err;

      gsl_fit_linear_est (xf, c0, c1, cov00, cov01, cov11, &yf, &yf_err);

      printf("fit: %g %g\n", xf, yf);
      printf("hi : %g %g\n", xf, yf + yf_err);
      printf("lo : %g %g\n", xf, yf - yf_err);
    }


  if (ok)
    return 0;
  return -1;
}