Re: [Help-gsl] Where oh where did the correlation cofficient go?

[Jack Denman]

At 02:25 AM 5/23/2007, you wrote:
>At Mon, 21 May 2007 23:11:15 -0700,
>Jack Denman wrote:
>> I discovered this great library, but in the linear regression code I
>> cannot find a way to get the correlation coefficient "rho" or r that
>> indicates the degree if fit. Am I missing something?
>
>Hello,
>
>We return the covariance matrix which gives more
>information, you can calculate the correlation
>coefficient from it as rho = c11/sqrt(c00*c11).
>
>-- 
>Brian Gough

I must have done something wrong. I created the results of an ideal linear 
function of y = mx +b and did the regression on the values and calculated the 
value of rho. I did this to check the accuracy of my code. The answer was off 
by a factor of roughly fifty ( rho = 0.02 ). I used both unweighted and 
neutrally weighted function calls. Below is listed the results and the code.

Where did I go wrong?

BTW: Thanks for the book advertisement.

Jack Denman
Fullerton California

*********************RESULTS************************************************************

# best fit: Y = 8.2 + 0.45 X
# covariance matrix:
# [ 0.0490741, -0.000925926  -0.000925926, 2.34412e-05]
# chisq = 1.54727e-28
r = 0.0218556
# best fit: Y = 8.2 + 0.45 X
# covariance matrix:
# [ 9.35661e-32, -1.7654e-33  -1.7654e-33, 4.46936e-35]
# sumsq = 1.48717e-28
r = 0.0218556
fit: -23.7 -2.465
hi : -23.7 -2.465
lo : -23.7 -2.465
fit: -22.91 -2.1095
hi : -22.91 -2.1095
lo : -22.91 -2.1095
fit: -22.12 -1.754
hi : -22.12 -1.754
lo : -22.12 -1.754
fit: -21.33 -1.3985
hi : -21.33 -1.3985
lo : -21.33 -1.3985
fit: -20.54 -1.043
hi : -20.54 -1.043
lo : -20.54 -1.043

*****************************CODE 
*********************************************************
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_fit.h>
/*
int gsl_fit_linear ( 
        const double *x, 
        const size_t xstride,
        const double *y,
        const size_t ystride,
        const size_t n,
        double *c0,
        double *c1,
        double *cov_00,
        double *cov_01,
        double *cov_11,
*/
int main ( void )
{
    int     i, n = 80;
    double  x[80] = {\
        0.000000,  1.000000,   2.000000, 3.000000,   4.000000,  5.000000,  
6.000000,  7.000000,  8.000000,  9.000000, \
        10.000000, 11.000000, 12.000000, 13.000000, 14.000000, 15.000000, 
16.000000, 17.000000, 18.000000, 19.000000, 20.000000,\
        21.000000, 22.000000, 23.000000, 24.000000, 25.000000, 26.000000, 
27.000000, 28.000000, 29.000000, 30.000000, 31.000000,\
        32.000000, 33.000000, 34.000000, 35.000000, 36.000000, 37.000000, 
38.000000, 39.000000, 40.000000, 41.000000, 42.000000,\
        43.000000, 44.000000, 45.000000, 46.000000, 47.000000, 48.000000, 
49.000000, 50.000000, 51.000000, 52.000000, 53.000000,\
        54.000000, 55.000000, 56.000000, 57.000000, 58.000000, 59.000000, 
60.000000, 61.000000, 62.000000, 63.000000, 64.000000,\
        65.000000, 66.000000, 67.000000, 68.000000, 69.000000, 70.000000, 
71.000000, 72.000000, 73.000000, 74.000000, 75.000000,\
        76.000000, 77.000000, 78.000000, 79.000000};
    double y [80]= {\
         8.200000,  8.650000,  9.100000,  9.550000, 10.000000, 10.450000, 
10.900000, 11.350000, 11.800000, 12.250000, 12.700000, 13.150000,\
        13.600000, 14.050000, 14.500000, 14.950000, 15.400000, 15.850000, 
16.300000, 16.750000, 17.200000, 17.650000, 18.100000, 18.550000,\
        19.000000, 19.450000, 19.900000, 20.350000, 20.800000, 21.250000, 
21.700000, 22.150000, 22.600000, 23.050000, 23.500000, 23.950000,\
        24.400000, 24.850000, 25.300000, 25.750000, 26.200000, 26.650000, 
27.100000, 27.550000, 28.000000, 28.450000, 28.900000, 29.350000,\
        29.800000, 30.250000, 30.700000, 31.150000, 31.600000, 32.050000, 
32.500000, 32.950000, 33.400000, 33.850000, 34.300000, 34.750000,\
        35.200000, 35.650000, 36.100000, 36.550000, 37.000000, 37.450000, 
37.900000, 38.350000, 38.800000, 39.250000, 39.700000, 40.150000,\
        40.600000, 41.050000, 41.500000, 41.950000, 42.400000, 42.850000, 
43.300000, 43.750000}; 
        
    double  w[80];
    for  ( i = 0; i < 80; i++ )
        w[i] = 1.0;
        
        
        
#if 0
    int     i, n = 4;
    double  x[4] = { 1970, 1980, 1990, 2000 };
    double  y[4] = { 12, 11, 14, 13 };
#endif
    double  c0, c1, cov00, cov01, cov11, chisq;
    double  cov_00, cov_01, cov_11, sumsq;
    size_t  xstride, ystride;
    int rc;
    gsl_fit_wlinear ( x, 1, w, 1, y, 1, n, &c0, &c1, &cov00, &cov01, &cov11, 
&chisq );
    printf ( "# best fit: Y = %g + %g X\n", c0, c1 );
    printf ( "# covariance matrix:\n" );
    printf ( "# [ %g, %g  %g, %g]\n", cov00, cov01, cov01, cov11 );
    printf ( "# chisq = %g\n", chisq );
    printf ("r = %g\n", cov11 / sqrt (cov00 * cov11 ) );
    //for ( i = 0; i < n; i++ )
//      printf ( "data: %g %g %g\n", x[i], y[i], 1/sqrt(w[i]) );
    gsl_fit_linear  ( x, 1, y, 1, n, &c0, &c1, &cov_00, &cov_01, &cov_11, 
&sumsq );
    printf ( "# best fit: Y = %g + %g X\n", c0, c1 );
    printf ( "# covariance matrix:\n" );
    printf ( "# [ %g, %g  %g, %g]\n", cov_00, cov_01, cov_01, cov_11 );
    printf ( "# sumsq = %g\n", sumsq );
    printf ("r = %g\n", cov_11 / sqrt (cov_00 * cov_11 ) );
    
    
//    for ( i = 0; i < n; i++ )
//      printf ( "data: %g %g\n", x[i], y[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 ); */
        gsl_fit_linear_est ( xf, c0, c1, cov_00, cov_01, cov_11, &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 );
    }
    return 0;
}