Function lfit in numerical recipes, providing a test function
- by Simon Walker
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 coe?cients 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
