#### 39.12.1 Exponential Fitting Example

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;
    }