Non-linear least squares fitting: Levemberg-Marquardts
Given a group of M pairs of experimental data (x1, y1), (x2, y2), ... , (xM, yM), having to fulfill certain function, and which depend on the variable x as well as on a group of N parameters aj
If the dependency of the aj parameters is non-linear, the y(x) function is expanded by means of a Taylor series of first order around the point aj0 (the initial values of the coefficients)
where
and yj represents the value of the function for aj.
The least squares method requires to minimize the expression:
with respect to each one of the coefficients, i.e.:
which leads to the following system of equations:
which can be written in matrix notation:
where:
The method of Marquardt exchanges matrix α in the previous equation by matrix α':
The files for this calculation are in folder breitwigner/3_Marquardt
. It uses the Numerical Recipes libraries mentioned above, plus covsrt.c
and gausj.c
.
The program marquardt.c
reads the experimental data from data.txt
and writes the results in an output.txt
file. The format of this file is given in three columns, where the first represents the energy, the second, the cross-section measured at that energy, and the third, the uncertainties. The values are taken from the experimental data of neutron scattering cross-sections provided in the main page.
In this case we have 3 parameters, σ, Er and γ. We have to provide initial values for them to start the iteration, and they will be adjusted by the mrqmin()
function. The program fits all the data points by default, with ia[i] = 1
.
/***************************************************
* Initial values for the parameter vectors
* For example:
* a[1] = 70000.0;
* a[2] = 75.0;
* a[3] = 60.0;
***************************************************/
for (i=1; i<=N; i++) {
printf("Enter initial value of parameter %d: ", i);
scanf("%lf", &a[i]);
}
// Fit all
for (i=1; i<=N; i++)
ia[i] = 1;
Once we have all the inputs, we can call the mrqmin()
function (you may want to check the definition again).
/***************************************************
* ITERATION LOOP
***************************************************/
char answer;
do {
printf("\nHow many iterations?: ");
scanf("%d", &k);
for(i=1; i<=k; i++) {
iter++;
ochisq = chisq;
mrqmin(x, y, sig, M, a, ia, N, covar, alpha, &chisq, funcs, &alamda);
if (chisq > ochisq)
n=0;
else if (fabs(ochisq - chisq) < 0.1)
n++;
if(n >= 4)
break;
}
if(n < 4) {
printf("\nIterate again? (Y/N)?: ");
scanf("%c", &answer);
}
} while (answer == 's' || answer == 'S');
Finally, we can use the user defined function func()
to recalculate the experimental points (see func()
definition ).
printf( "\nData points recalculated from the fit:\n");
fprintf(fp, "\nData points recalculated from the fit:\n");
for (i=1; i<=M; i++) {
(*funcs)(x[i], a, &ymod[i], dyda, N);
sum = 0.0;
for (j=1; j<=N; j++)
sum += pow(dyda[j]*sqrt(covar[j][j]), 2);
printf( "\nSigma[%d] = %3.3lf\t\tdSigma[%d] = %3.3lf", i, ymod[i], i, sqrt(sum));
fprintf(fp, "\nSigma[%d] = %3.3lf\t\tdSigma[%d] = %3.3lf", i, ymod[i], i, sqrt(sum));
}
Results
The following values with their errors are obtained for the parameters of interest:
- σ0 = 66900.3
- Er = 77.5 MeV
- γ = 56.2 MeV