Greetings!
I have been trying to use gsl for some linear least squares
fitting, but for some reason I'm getting bad results.
Instead of returning the results of the fit, gsl_fit_linear()
returns -nan for all parameters. I'm not sure if I'm passing something
wrong to it, but I do know that there is no problem with the data (I've
ran it through other LS fitting procedures and results are consistent).
Attached are the source and an example data file. I've compiled
and ran with:
$ gcc -I/usr/include/gsl/ -lm -lgsl -lgslcblas -o teste_bet teste_bet.c
$ ./teste_bet data.dat
--
--
[]s,
Francisco M Neto-3.14699817384996E-05 4.09163756453624
4.07463535162317E-05 8.09641771055868
8.43264717798022E-05 28.2838390868175
0.00018798259043 66.8821177184521
0.00030092780419 71.881475368167
0.00046269415648 77.1289379089803
0.00055319581577 82.3067049142331
0.00073340490794 84.0925925488574
0.00079172642877 87.1782054368395
0.00106511299186 89.6083469704575
0.0012225076642 93.4451737206919
0.00160956002107 96.5966938587599
0.0018399821398 100.524550335042
0.0023146733882 104.38666955937
0.00300479604127 108.671730117131
0.00357303476627 112.400630270797
0.00444696982072 117.227524747886
0.00576705677146 122.004429070564
0.00718977438218 126.808396744222
0.0089289235829 131.955525144443
0.01077396805166 135.931230923791
0.01338169911633 141.153959875754
0.01760289341922 147.643170077296
0.02399056709581 155.880429521727
0.03185636554344 163.612166329895
0.03761380968281 168.318128484672
0.04504245151261 173.469057637336
0.05160420469665 177.857318295557
0.05700608358335 181.125025882971
0.06573572201811 185.879315463347
0.07254411320494 189.298863447612
0.07943660024635 192.470406857215
0.08630575663586 195.348252051481
0.0929089010198 198.225735225978
0.09976326687852 201.088784192443
0.10667499722023 203.787515758126
0.11948249650093 208.527244894866
0.14304333823935 216.508971798891
0.16915811680024 224.593361511492
0.19403117331301 231.7222283824
0.21883175193984 238.725516919965
0.24413051941145 245.611693523167
0.26864470363083 252.083392881103
0.29387449819356 258.575751712931
0.31791810413123 264.77314373454
0.34218982622314 271.051892151926
0.36220734105375 276.089002880462
0.38588320561528 282.322676586783
0.4106427694869 289.033307635841
0.43484981744396 295.729803052409
0.45893371052068 302.499665377277
0.48315166150882 309.543152413904
0.50755746316707 316.885947661475
0.53171349985555 324.48069434984
0.55599835203823 332.452951747861
0.57918790488814 340.57540762982
0.60360059244569 349.592858469625
0.62765322297741 359.01986812748
0.65135097032746 369.418462996741
0.67626305266922 381.373445868837
0.7000888588966 394.225604825983
0.70617784085263 397.800803104449
0.71025246360704 400.306112219663
0.71489684133586 403.184253054187
0.71948763798889 406.140288418243
0.72396491285022 409.215888535516
0.72863727878662 412.713008977103
0.7332194809001 416.675363782455
0.73757485667776 421.861767927571
0.74192908408716 430.362362945718
0.74601430465305 444.420400921572
0.75000391682773 467.273482183065
0.75379164961459 494.833343070218
0.75777196653392 525.080288294786
0.76192341456643 558.275409865907
0.7668598805338 592.70306599444
0.77114448001604 627.811102273906
0.77557038853593 664.048712652327
0.781590282284 699.964187719685
0.80119953417553 718.759907300559
0.8296574674607 726.809353933258
0.85062236435354 732.178729396104
0.85393986715313 733.027199440747
0.87334262633399 737.019872910495
0.88968204475573 739.830163214467
0.90739318296275 742.994989119476
0.92614263047979 746.514948649767
0.94381537891845 750.141396358427
0.96130737847078 754.725900782116
0.9710888622052 758.659843335282
0.98083226283641 764.492162151857
0.98840131246824 773.889616411278
0.99422319914231 789.957479214035
0.98668865285944 774.100165526189
0.96864595705195 759.257813041533
0.94349447124659 751.328844541608
0.92242503488606 747.061018954581
0.90430611085806 743.711650970592
0.88694167540097 740.79010672564
0.87008428658249 737.709095457114
0.86863771258244 737.430591839287
0.85178747377444 734.071830545185
0.83530449736255 730.311279502191
0.81701602679231 726.51010517774
0.80097712396177 723.406608863467
0.7842054845535 720.080673831493
0.76706672671769 716.45551789223
0.74974739857146 712.814285549007
0.74494088843019 711.610723821362
0.73972090708891 710.237463737819
0.73477367982084 708.916490797052
0.73007573615619 707.372139545577
0.72513293803212 705.002532244208
0.72001312912967 701.252431798027
0.71577963020039 695.801476584886
0.71223244921034 685.86259233302
0.70862061626817 668.027545315774
0.7056795123899 638.423844763533
0.69539440098059 522.935698593525
0.67616712908857 402.316876348799
0.64816886949234 369.358013213783
0.62177224350371 357.896702974209
0.59874571159119 348.569168753684
0.58779135618897 344.517551492138
0.56738879130956 337.284464107615
0.54662609402197 330.287385281973
0.52647285898516 323.63934466543
0.50621165774291 317.097270013
0.48547062954318 310.808601008101
0.46552057111792 305.033642518225
0.44482061361254 298.914577272251
0.42365124029891 292.938493013621
0.40380327341393 287.510158953291
0.38333334360807 282.059852533332
0.36333678390995 276.810917398316
0.34261156088373 271.235566663875
0.32155042986013 265.862056983182
0.30181886000268 260.824256526279
0.28076844876272 255.417363041364
0.26043701451108 250.134817727256
0.23989397781208 244.709629631459
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl_fit.h>
/* Maximum name size + 1. */
#define MAX_NAME_SZ 256
#define MAX_LINE_SIZE 256
int file_lines (FILE *file);
int file_read (FILE *file, int n);
int main (int argc, char *argv[]) {
int lines, i;
FILE *input;
char *name, line[MAX_LINE_SIZE];
double **data, **bet;
double a, b, cov00, cov01, cov11, chisq; // Fit variables
int n, status;
if (!(name = malloc(sizeof(char)*MAX_NAME_SZ))) {
printf("No memory!\n");
return 1;
}
if (argc == 1) { // No extra arguments on the command line
printf("Name of the input file: ");
fgets(name, MAX_NAME_SZ, stdin);
} else if (argc == 2) {
strcpy(name, argv[1]);
}
else {
printf ("Too many arguments! Exiting.\n");
exit(1);
}
if ((strlen(name) > 0) && (name[strlen (name) - 1] == '\n'))
name[strlen (name) - 1] = '\0';
//printf ("Loading file: %s\n", name);
if (!(input=fopen(name, "r")))
exit(1);
lines = file_lines (input);
// printf ("Lines in the file: %d\n", lines);
if(!(data = malloc(2 * sizeof(double *))))
exit(1);
for (i=0; i<2; i++)
if(!(data[i] = malloc((lines) * sizeof(double))))
exit(1);
for (i=0; i<lines; i++) {
fgets(line, MAX_LINE_SIZE, input);
line[strlen (line) - 1] = '\0';
data[0][i] = atof(strtok(line, " \t"));
data[1][i] = atof(strtok(NULL, " \t"));
}
// File loaded into data[][].
if(!(bet = malloc(2 * sizeof(double *))))
exit(1);
for (i=0; i<=2; i++)
if(!(bet[i] = malloc((lines) * sizeof(double))))
exit(1);
n = 0;
for (i=0; i<lines; i++) {
if (data[0][i] >= 0.05 && data[0][i] <= 0.3) {
bet[0][i] = data[0][i];
bet[1][i] = 1/(data[1][i]*((1/data[0][i])-1)); // y' = 1/(y*(x-1))
bet[2][i] = 1/(bet[1][i]*0.05); // weight for the fits, based on 5% estimated error
printf ("%.8e\t%.8e\t%.8e\n", bet[0][i], bet[1][i], bet[2][i]);
n++;
}
}
// Linear fit to BET-transformed data from 0.05 < x < 0.3 (y' = a + bx')
printf ("Performing fit with %d points...", n);
status = gsl_fit_linear (bet[0], 1, bet[1], 1, n, &a, &b, &cov00, &cov01, &cov11, &chisq);
if (status) printf ("not ");
printf ("done! (status = %d)\n", status);
printf ("====== BET results ======\n");
printf ("y = %lf + %lf x\n", a, b);
printf ("[%g, %g\n %g, %g]\n",
cov00, cov01, cov01, cov11);
printf ("Chi Square = %lf\n", chisq);
printf ("=========================\n");
fclose (input);
free (input);
free (name);
free (data);
free (bet);
return 0;
}
int file_lines (FILE *file) { // Counts the lines in a text file
int i=0, c;
while (!feof(file)) {
c = fgetc(file);
if (c == '\n')
i++;
}
rewind(file);
return i;
}
