[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Polynomial Regression

From: Jack Denman
Subject: [Help-gsl] Polynomial Regression
Date: Wed, 20 Jun 2007 09:34:29 -0700

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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]