title: Fit a gaussian with GSL GNU (levmar algo)
body: <p>I'm using Nonlinear Least-Squares Fitting og GSL library to fit a gaussian.
  I used the example provided from <a href="https://www.gnu.org/software/gsl/doc/html/nls.html" rel="nofollow noreferrer">https://www.gnu.org/software/gsl/doc/html/nls.html</a> (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</p>
  <p>So I made the following modifications :
  I added the offset in the gaussian fct:</p>
  <pre><code>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);
  }
  </code></pre>
  <p>and I updated the jacobian func</p>
  <pre><code>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 &lt; d-&gt;n; ++i)
      {
        double ti = d-&gt;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;
  }
  </code></pre>
  <p>When I run an example, the algo is not able to fit the data (it perfroms one iteration):</p>
  <p>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</p>
  <p>NITER         = 1</p>
  <p>NFEV          = 33</p>
  <p>NJEV          = 1</p>
  <p>NAEV          = 0</p>
  <p>initial cost  = 1.006561493040e+002</p>
  <p>final cost    = 1.006561493040e+002</p>
  <p>final x       = (1.000000000000e+000, 0.000000000000e+000, 1.000000e+000, 0.000000e+000)</p>
  <p>final cond(J) = 9.099099513949e+001</p>
  <p>reason for stopping: small gradient</p>
  <p>Here the values of the example:</p>
  <pre><code>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);
  </code></pre>
  <p>Any idea?
  I can provide the full example</p>