Next: Nonlinear Least-Squares Geodesic Acceleration Example, Up: Nonlinear Least-Squares Examples [Index]
The following example program fits a weighted exponential model with
background to experimental data, Y = A \exp(-\lambda t) + b. The
first part of the program sets up the functions expb_f and
expb_df to calculate the model and its Jacobian. The appropriate
fitting function is given by,
f_i = (A \exp(-\lambda t_i) + b) - y_i
where we have chosen t_i = i. The Jacobian matrix J is the derivative of these functions with respect to the three parameters (A, \lambda, b). It is given by,
J_{ij} = d f_i / d x_j
where x_0 = A, x_1 = \lambda and x_2 = b. The ith row of the Jacobian is therefore
The main part of the program sets up a Levenberg-Marquardt solver and some simulated random data. The data uses the known parameters (5.0,0.1,1.0) combined with Gaussian noise (standard deviation = 0.1) over a range of 40 timesteps. The initial guess for the parameters is chosen as (1.0, 1.0, 0.0). The iteration terminates when the relative change in x is smaller than 10^{-8}, or when the magnitude of the gradient falls below 10^{-8}. Here are the results of running the program:
iter 0: A = 1.0000, lambda = 1.0000, b = 0.0000, cond(J) = inf, |f(x)| = 62.2029 iter 1: A = 1.2196, lambda = 0.3663, b = 0.0436, cond(J) = 53.6368, |f(x)| = 59.8062 iter 2: A = 1.6062, lambda = 0.1506, b = 0.1054, cond(J) = 23.8178, |f(x)| = 53.9039 iter 3: A = 2.4528, lambda = 0.0583, b = 0.2470, cond(J) = 20.0493, |f(x)| = 28.8039 iter 4: A = 2.9723, lambda = 0.0494, b = 0.3727, cond(J) = 94.5601, |f(x)| = 15.3252 iter 5: A = 3.3473, lambda = 0.0477, b = 0.4410, cond(J) = 229.3627, |f(x)| = 10.7511 iter 6: A = 3.6690, lambda = 0.0508, b = 0.4617, cond(J) = 298.3589, |f(x)| = 9.7373 iter 7: A = 3.9907, lambda = 0.0580, b = 0.5433, cond(J) = 250.0194, |f(x)| = 8.7661 iter 8: A = 4.2353, lambda = 0.0731, b = 0.7989, cond(J) = 154.8571, |f(x)| = 7.4299 iter 9: A = 4.6573, lambda = 0.0958, b = 1.0302, cond(J) = 140.2265, |f(x)| = 6.1893 iter 10: A = 5.0138, lambda = 0.1060, b = 1.0329, cond(J) = 109.4141, |f(x)| = 5.4961 iter 11: A = 5.1505, lambda = 0.1103, b = 1.0497, cond(J) = 100.8762, |f(x)| = 5.4552 iter 12: A = 5.1724, lambda = 0.1110, b = 1.0526, cond(J) = 97.3403, |f(x)| = 5.4542 iter 13: A = 5.1737, lambda = 0.1110, b = 1.0528, cond(J) = 96.7136, |f(x)| = 5.4542 iter 14: A = 5.1738, lambda = 0.1110, b = 1.0528, cond(J) = 96.6678, |f(x)| = 5.4542 iter 15: A = 5.1738, lambda = 0.1110, b = 1.0528, cond(J) = 96.6663, |f(x)| = 5.4542 iter 16: A = 5.1738, lambda = 0.1110, b = 1.0528, cond(J) = 96.6663, |f(x)| = 5.4542 summary from method 'trust-region/levenberg-marquardt' number of iterations: 16 function evaluations: 23 Jacobian evaluations: 17 reason for stopping: small step size initial |f(x)| = 62.202928 final |f(x)| = 5.454180 chisq/dof = 0.804002 A = 5.17379 +/- 0.27938 lambda = 0.11104 +/- 0.00817 b = 1.05283 +/- 0.05365 status = success
The approximate values of the parameters are found correctly, and the chi-squared value indicates a good fit (the chi-squared per degree of freedom is approximately 1). In this case the errors on the parameters can be estimated from the square roots of the diagonal elements of the covariance matrix. If the chi-squared value shows a poor fit (i.e. chi^2/dof >> 1) then the error estimates obtained from the covariance matrix will be too small. In the example program the error estimates are multiplied by \sqrt{\chi^2/dof} in this case, a common way of increasing the errors for a poor fit. Note that a poor fit will result from the use of an inappropriate model, and the scaled error estimates may then be outside the range of validity for Gaussian errors.
Additionally, we see that the condition number of J(x) stays reasonably small throughout the iteration. This indicates we could safely switch to the Cholesky solver for speed improvement, although this particular system is too small to really benefit.
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>
/* number of data points to fit */
#define N 40
struct data {
size_t n;
double * y;
};
int
expb_f (const gsl_vector * x, void *data,
gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *y = ((struct data *)data)->y;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
double b = gsl_vector_get (x, 2);
size_t i;
for (i = 0; i < n; i++)
{
/* Model Yi = A * exp(-lambda * i) + b */
double t = i;
double Yi = A * exp (-lambda * t) + b;
gsl_vector_set (f, i, Yi - y[i]);
}
return GSL_SUCCESS;
}
int
expb_df (const gsl_vector * x, void *data,
gsl_matrix * J)
{
size_t n = ((struct data *)data)->n;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
size_t i;
for (i = 0; i < n; i++)
{
/* Jacobian matrix J(i,j) = dfi / dxj, */
/* where fi = (Yi - yi)/sigma[i], */
/* Yi = A * exp(-lambda * i) + b */
/* and the xj are the parameters (A,lambda,b) */
double t = i;
double e = exp(-lambda * t);
gsl_matrix_set (J, i, 0, e);
gsl_matrix_set (J, i, 1, -t * A * e);
gsl_matrix_set (J, i, 2, 1.0);
}
return GSL_SUCCESS;
}
void
callback(const size_t iter, void *params,
const gsl_multifit_nlinear_workspace *w)
{
gsl_vector *f = gsl_multifit_nlinear_residual(w);
gsl_vector *x = gsl_multifit_nlinear_position(w);
double rcond;
/* compute reciprocal condition number of J(x) */
gsl_multifit_nlinear_rcond(&rcond, w);
fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n",
iter,
gsl_vector_get(x, 0),
gsl_vector_get(x, 1),
gsl_vector_get(x, 2),
1.0 / rcond,
gsl_blas_dnrm2(f));
}
int
main (void)
{
const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
gsl_multifit_nlinear_workspace *w;
gsl_multifit_nlinear_fdf fdf;
gsl_multifit_nlinear_parameters fdf_params =
gsl_multifit_nlinear_default_parameters();
const size_t n = N;
const size_t p = 3;
gsl_vector *f;
gsl_matrix *J;
gsl_matrix *covar = gsl_matrix_alloc (p, p);
double y[N], weights[N];
struct data d = { n, y };
double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
gsl_vector_view x = gsl_vector_view_array (x_init, p);
gsl_vector_view wts = gsl_vector_view_array(weights, n);
gsl_rng * r;
double chisq, chisq0;
int status, info;
size_t i;
const double xtol = 1e-8;
const double gtol = 1e-8;
const double ftol = 0.0;
gsl_rng_env_setup();
r = gsl_rng_alloc(gsl_rng_default);
/* define the function to be minimized */
fdf.f = expb_f;
fdf.df = expb_df; /* set to NULL for finite-difference Jacobian */
fdf.fvv = NULL; /* not using geodesic acceleration */
fdf.n = n;
fdf.p = p;
fdf.params = &d;
/* this is the data to be fitted */
for (i = 0; i < n; i++)
{
double t = i;
double yi = 1.0 + 5 * exp (-0.1 * t);
double si = 0.1 * yi;
double dy = gsl_ran_gaussian(r, si);
weights[i] = 1.0 / (si * si);
y[i] = yi + dy;
printf ("data: %zu %g %g\n", i, y[i], si);
};
/* allocate workspace with default parameters */
w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);
/* initialize solver with starting point and weights */
gsl_multifit_nlinear_winit (&x.vector, &wts.vector, &fdf, w);
/* compute initial cost function */
f = gsl_multifit_nlinear_residual(w);
gsl_blas_ddot(f, f, &chisq0);
/* solve the system with a maximum of 20 iterations */
status = gsl_multifit_nlinear_driver(20, xtol, gtol, ftol,
callback, NULL, &info, w);
/* compute covariance of best fit parameters */
J = gsl_multifit_nlinear_jac(w);
gsl_multifit_nlinear_covar (J, 0.0, covar);
/* compute final cost */
gsl_blas_ddot(f, f, &chisq);
#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
fprintf(stderr, "summary from method '%s/%s'\n",
gsl_multifit_nlinear_name(w),
gsl_multifit_nlinear_trs_name(w));
fprintf(stderr, "number of iterations: %zu\n",
gsl_multifit_nlinear_niter(w));
fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
fprintf(stderr, "Jacobian evaluations: %zu\n", fdf.nevaldf);
fprintf(stderr, "reason for stopping: %s\n",
(info == 1) ? "small step size" : "small gradient");
fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
fprintf(stderr, "final |f(x)| = %f\n", sqrt(chisq));
{
double dof = n - p;
double c = GSL_MAX_DBL(1, sqrt(chisq / dof));
fprintf(stderr, "chisq/dof = %g\n", chisq / dof);
fprintf (stderr, "A = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
fprintf (stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
fprintf (stderr, "b = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
}
fprintf (stderr, "status = %s\n", gsl_strerror (status));
gsl_multifit_nlinear_free (w);
gsl_matrix_free (covar);
gsl_rng_free (r);
return 0;
}