[Help-gsl] Polynomial Regression

[Jack Denman]

I calculated the linear regression on some financial data of the income
and finances of the association I belong to - I am a director. Some of
the correlation coefficients were poor so I made another attempt at a
2nd degree polynomial and received a significantly lower sum of the
squares (although with no correlation value) and a better fit just
observed from the data although I have not actually plotted the data
yet. I remember from the "bad old days" a FORTRAN package that called
statpack or statpac than started with linear regression and proceeded to
increase the degree of the polynomial until there was no further
reduction in the sum of the squares. Is there a way to run a cubic
regression. My second degree regression was with the following calls
from examples in the gsl HTML page.

Jack Denman

Note yy = array & xx = array of n count.
    mean = gsl_stats_mean (yy, 1,  n);
    std  = gsl_stats_sd_m (yy, 1, n, mean);

    X = gsl_matrix_alloc (n, 3);
    y = gsl_vector_alloc (n);
    w = gsl_vector_alloc (n);

    c = gsl_vector_alloc (3);
    cov = gsl_matrix_alloc (3, 3);

    ei = std;


  for (i = 0; i < n; i++)
    {

        xi = xx[i];
        yi = yy[i];

        printf ("%g %g +/- %g\n", xi, yi, ei);

        gsl_matrix_set (X, i, 0, 1.0);
        gsl_matrix_set (X, i, 1, xi);
        gsl_matrix_set (X, i, 2, xi * xi);

        gsl_vector_set (y, i, yi);
        gsl_vector_set (w, i, 1.0 / (ei * ei));
    }

    {
        gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc
(n, 3);
        gsl_multifit_wlinear (X, w, y, c, cov, &chisq, work);
        gsl_multifit_linear_free (work);
    }

#define C(i) (gsl_vector_get(c,(i)))
#define COV(i,j) (gsl_matrix_get(cov,(i),(j)))

    {
        printf ("# best fit: Y = %g + %g X + %g X^2\n", C (0), C (1), C
(2));
        printf ("# covariance matrix:\n");
        printf ("[ %+.5e, %+.5e, %+.5e  \n", COV (0, 0), COV (0, 1), COV
(0, 2));
        printf ("  %+.5e, %+.5e, %+.5e  \n", COV (1, 0), COV (1, 1), COV
(1, 2));
        printf ("  %+.5e, %+.5e, %+.5e ]\n", COV (2, 0), COV (2, 1), COV
(2, 2));
        printf ("# chisq = %g\n", chisq);
   }