Function lfit in numerical recipes, providing a test function

Question

Hi

I am trying to fit collected data to a polynomial equation and I found the lfit function from Numerical Recipes. I only have access to the second edition, so am using that.

I have read about the lfit function and its parameters, one of which is a function pointer, given in the documentation as

void (*funcs)(float, float [], int))

with the help

The user supplies a routine funcs(x,afunc,ma) that returns the ma basis functions evaluated at x = x in the array afunc[1..ma].

I am struggling to understand how this lfit function works. An example function I found is given below:

void fpoly(float x, float p[], int np)
/*Fitting routine for a polynomial of degree np-1, with coefficients in the array p[1..np].*/
{
    int j;
    p[1]=1.0;
    for (j=2;j<=np;j++)
        p[j]=p[j-1]*x;
}

When I run through the source code for the lfit function in gdb I can see no reference to the funcs pointer. When I try and fit a simple data set with the function, I get the following error message.

Numerical Recipes run-time error...
gaussj: Singular Matrix
...now exiting to system...

Clearly somehow a matrix is getting defined with all zeroes.

I am going to involve this function fitting in a large loop so using another language is not really an option. Hence why I am planning on using C/C++.

For reference, the test program is given here:

int main()
{

    float x[5] = {0., 0., 1., 2., 3.};
    float y[5] = {0., 0., 1.2, 3.9, 7.5};
    float sig[5] = {1., 1., 1., 1., 1.};

    int ndat = 4;
    int ma = 4; /* parameters in equation */
    float a[5] = {1, 1, 1, 0.1, 1.5};
    int ia[5] = {1, 1, 1, 1, 1};
    float **covar = matrix(1, ma, 1, ma);

    float chisq = 0;

    lfit(x,y,sig,ndat,a,ia,ma,covar,&chisq,fpoly);

    printf("%f\n", chisq);  
    free_matrix(covar, 1, ma, 1, ma);


    return 0;
}

Also confusing the issue, all the Numerical Recipes functions are 1 array-indexed so if anyone has corrections to my array declarations let me know also!

Cheers

edit: Ok just discovered the problem to this.

When copying the numerical recipes code, I accidentally commented out some real code thinking it was just a comment. I no longer get the singular value error

Answer 1

I don't know much about the lfit function but could you please add some fprintf to stderr at the beginning of the fpoly function and another at the end of the function to see where it actually crashes or if the parameters received are correct/or valid?

Sorry if this is not helpful but I assume that lfit function is working well and the problem is in the fpoly function.

Answer 2

I have since switched to the 3rd edition of NR which has a more understandable routine.