Francisco,
You have some errors in the way you handle your 2D arrays. I have
made a small modification to define new 1D arrays (x and y) and store
your fit data into those arrays. The program output now says:
Performing fit with 20 points...done! (status = 0)
====== BET results ======
y = 0.000022 + 0.005321 x
[2.54531e-11, -1.26944e-10
-1.26944e-10, 8.07463e-10]
Chi Square = 0.000000
=========================
My changes are attached - you can find the lines I changed with the "PA"
keyword
On 05/03/2018 04:13 AM, Francisco M Neto wrote:
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
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/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 *x, *y; // PA
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);
// PA: allocate x and y
x = malloc(lines * sizeof(double));
y = malloc(lines * sizeof(double));
n = 0;
for (i=0; i<lines; i++) {
if (data[0][i] >= 0.05 && data[0][i] <= 0.3) {
// PA: set x and y vectors
x[n] = data[0][i];
y[n] = 1/(data[1][i]*((1/data[0][i])-1));
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')
#if 1
{
printf ("Performing fit with %d points...", n);
status = gsl_fit_linear (x, 1, y, 1, n, &a, &b, &cov00, &cov01, &cov11, &chisq);
if (status) printf ("not ");
printf ("done! (status = %d)\n", status);
}
#else
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);
#endif
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 (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;
}
