least squares fitting in c

isamuede · 09-05-2009, 08:27 PM

Hi there, im using the gsl library and try to play with the least squares fitting routines so I test the next code but when compiling it it shows me warnings and at the end it does not show me the y=mx+b thing
only shows a nan string. Any help qould be apreciated thanks!!!

Code:

#include <stdio.h>
#include <gsl/gsl_fit.h>
int
main (void)
{
  int i, n = 13;
  double x[5] = { 1.08, 1.14, 1.2, 1.26, 1.3, 1.36, 1.42, 1.5, 1.6, 1.7, 1.84, 2.02, 2.24 };
  double y[5] = {   0.110, 0.105, 0.100, 0.095, 0.090, 0.085, 0.080, 0.075, 0.070, 0.065, 0.060, 0.055, 0.050 };
  double c0, c1, cov00, cov01, cov11, sumsq;
  gsl_fit_linear (x, 1, y, 1, n,
                   &c0, &c1, &cov00, &cov01, &cov11,
                   &sumsq);
  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 ("# sumsq = %g\n", sumsq);
  for (i = 0; i < n; i++)
    printf ("data: %g %g %g\n",
                   x[i], y[i], 1/sqrt(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);
          }
       return 0;
     }

when using
gcc -Wall -I/usr/local/include -c thomson.c
it says
thomson.c:8: warning: (near initialization for ‘y’)
and then
gcc -L/usr/local/lib thomson.o -lgsl -lgslcblas -lm

and running it
./a.out
# best fit: Y = nan + nan X
# covariance matrix:
# [ nan, nan
# nan, nan]
# sumsq = nan
data: 1.08 0.11 inf
data: 1.14 0.105 1
data: 1.2 0.1 0.707107
data: 1.26 0.095 0.57735
data: 1.3 0.09 0.5
data: 2.75859e-313 1.08 0.447214
data: 1.23694e+13 1.14 0.408248
data: 6.64592e-316 1.2 0.377964
data: 4.70127e+12 1.26 0.353553
data: -0.0254148 1.3 0.333333
data: -9.35615e-39 2.75859e-313 0.316228
data: nan 1.23694e+13 0.301511
data: 4.8537e-270 6.64592e-316 0.288675

fit: 1.404 nan
hi : 1.404 nan
lo : 1.404 nan
fit: 1.3932 nan
hi : 1.3932 nan
lo : 1.3932 nan
fit: 1.3824 nan
hi : 1.3824 nan
lo : 1.3824 nan
fit: 1.3716 nan
hi : 1.3716 nan
lo : 1.3716 nan
fit: 1.3608 nan
hi : 1.3608 nan
lo : 1.3608 nan
fit: 1.35 nan
hi : 1.35 nan
lo : 1.35 nan
fit: 1.3392 nan
hi : 1.3392 nan

thanks again

paulsm4 · 09-05-2009, 09:48 PM

Hi -

The first problem I see is this:

Code:

/* Bad: declare static array of "5" ... then initialize with "13"! */
  double x[5] = { 1.08, 1.14, 1.2, 1.26, 1.3, 1.36, 1.42, 1.5, 1.6, 1.7, 1.84, 2.02, 2.24 };
  double y[5] = {   0.110, 0.105, 0.100, 0.095, 0.090, 0.085, 0.080, 0.075, 0.070, 0.065, 0.060, 0.055, 0.050 };

Code:

/* Better: static array, no explicit size; initialize with as many/few as you want */
  double x[] = {  1.08,  1.14,   1.2,  1.26,   1.3,  1.36,  1.42,   1.5,   1.6,   1.7,  1.84,  2.02,  2.24 };
  double y[] = { 0.110, 0.105, 0.100, 0.095, 0.090, 0.085, 0.080, 0.075, 0.070, 0.065, 0.060, 0.055, 0.050 };

Code:

/* Better still ... */
#include <stdio.h>
#include <gsl/gsl_fit.h>

#define NELMS(A) (sizeof(A) / sizeof(A[0]) )

int
main (void)
{
  double x[] = {  1.08,  1.14,  1.2,   1.26,   1.3,  1.36,  1.42,   1.5,   1.6,   1.7,  1.84,  2.02,  2.24 };
  double y[] = { 0.110, 0.105, 0.100, 0.095, 0.090, 0.085, 0.080, 0.075, 0.070, 0.065, 0.060, 0.055, 0.050 };
  int i = 0, n = NELMS(x);
  double c0 = 0, c1 = 0, cov00 = 0, cov01 = 0, cov11 = 0, sumsq = 0;
  ...

'Hope that helps .. PSM

isamuede · 09-05-2009, 10:08 PM

whoa!! You see the effects of the distraction , thanks paulsm4!!!!

isamuede · 09-05-2009, 10:10 PM

Code:

#define NELMS(A) (sizeof(A) / sizeof(A[0])

Very smart!!!!!
Thanks again

isamuede · 09-05-2009, 10:15 PM

Code:

#define NELMS(A) (sizeof(A) / sizeof(A[0])

Very smart!!!!!
Thanks again

Sergei Steshenko · 09-05-2009, 11:26 PM

Quote:

Originally Posted by isamuede

Code:

#define NELMS(A) (sizeof(A) / sizeof(A[0])

Very smart!!!!!
Thanks again

Yes, it's smart, and it's in "C" FAQ for ages.

wje_lq · 09-06-2009, 04:29 AM

Quote:

Originally Posted by Sergei Steshenko

Yes, it's smart, and it's in "C" FAQ for ages.

Exactly right.

But to be slightly more helpful, you can find it here.

Browse through the whole document, isamuede. It's likely to open your eyes.

isamuede · 09-06-2009, 11:13 AM

Thanks wje_lq , im checking it right now~~!!! You are very kind!!!