Fit a gaussian with GSL GNU (levmar algo)

Question

I'm using Nonlinear Least-Squares Fitting og GSL library to fit a gaussian. I used the example provided from https://www.gnu.org/software/gsl/doc/html/nls.html (Example 2). It fits succesfully a gaussian with 3 parameters (a, b, c) that represents resp the amplitude, mean, and width (even without the Geodesic Acceleration). Now I would like to add a parameter (an offset) as the function to fit is now: a * exp(-0.5 (t - b)² / c²) + d

So I made the following modifications : I added the offset in the gaussian fct:

gaussian(const double a, const double b, const double c, const double d, const double t)
{
  const double z = (t - b) / c;
  return (a * exp(-0.5 * z * z) + d);
}

and I updated the jacobian func

int func_df (const gsl_vector * x, void *params, gsl_matrix * J)
{
  struct data *d = (struct data *) params;
  double a = gsl_vector_get(x, 0);
  double b = gsl_vector_get(x, 1);
  double c = gsl_vector_get(x, 2);
  size_t i;

  for (i = 0; i < d->n; ++i)
    {
      double ti = d->t[i];
      double zi = (ti - b) / c;
      double ei = exp(-0.5 * zi * zi);

      gsl_matrix_set(J, i, 0, -ei);
      gsl_matrix_set(J, i, 1, -(a / c) * ei * zi);
      gsl_matrix_set(J, i, 2, -(a / c) * ei * zi * zi);
      gsl_matrix_set(J, i, 3, 1);
    }

  return GSL_SUCCESS;
}

When I run an example, the algo is not able to fit the data (it perfroms one iteration):

iter 0: a = 1.0000, b = 0.0000, c = 1.0000, d = 0.0000, |a|/|v| = 0.0000 cond(J) = 1.#INF, |f(x)| = 10.0328

NITER = 1

NFEV = 33

NJEV = 1

NAEV = 0

initial cost = 1.006561493040e+002

final cost = 1.006561493040e+002

final x = (1.000000000000e+000, 0.000000000000e+000, 1.000000e+000, 0.000000e+000)

final cond(J) = 9.099099513949e+001

reason for stopping: small gradient

Here the values of the example:

const size_t n = 100;  /* number of data points to fit */
const size_t p = 4;    /* number of model parameters */

/*Gaussian to be fitted*/
const double a = 5.0;  /* amplitude */
const double b = 0.4;  /* center */
const double c = 0.15; /* width */
const double dd = 1.0; /* offset */

/* starting point */
gsl_vector_set(x, 0, 1.0);
gsl_vector_set(x, 1, 0.0);
gsl_vector_set(x, 2, 1.0);
gsl_vector_set(x, 3, 0.0);

Any idea? I can provide the full example

Fit a gaussian with GSL GNU (levmar algo)

0 Answers