devernay / cminpack

A C/C++ rewrite of the MINPACK software (originally in FORTRAN) for solving nonlinear equations and nonlinear least squares problems
http://devernay.free.fr/hacks/cminpack/
145 stars 63 forks source link

[question] I don't get any result when using cminpack routines for the least squares nonlinear problem. #43

Closed supernova-III closed 3 years ago

supernova-III commented 3 years ago

Hi everyone! I want to use cmipack to fit my data to different models (least squares nonlinear problem) but actually I see that routine returns 1 and initial guess vector is not changed, but the fitting problem is quite simple. I think that probably I'm doing something wrong.

My fitting problem is the following. Find parameters of gaussian function (amplitude, mu, and sigma) to fit the data:

  // abscissa
  const std::vector<double> x1{
    -3.01969624096814,
    -2.20911406128415,
    -1.39853188160016,
    -0.587949701916172,
    0.222632477767818,
    1.03321465745181,
    1.8437968371358,
    2.65437901681979,
    3.46496119650378
  };

  // data
  const std::vector<double> y1{
    22, 183, 917, 2446, 3239, 2263, 778, 140, 12
  };

It's obvious gaussian function as we can see. In cmipack we can use several functions to fit the data. I've tried lmdif1, when the derivatives are calculated automatically using forward finite difference approximation. In this case I need just define the function to minimize:

struct FData
{
  std::vector<double> x;
  std::vector<double> dataToModel;
};
int gaussian_cminpack(void* p, int m, int n, const double* x,
                       double* fvec, int flag)
{
  auto myData = reinterpret_cast<FData*>(p);

  for (int i = 0; i < m; ++i)
  {
    double y = x[2]  * exp(-0.5 * sqr((myData->x[i] - x[0]) / x[1])) / (sqrt(2 * PI) * x[1]);
    fvec[i] = abs(y - sqr(myData->dataToModel[i]));
  }
  return 0;
}

Then I just following the documentation and

  double iGuess[]{ 0.01, 0.6, 3239 };

  FData data{ x1, y1 };
  double fvec[9];
  int iwa[3];
  double wa[51];
  int lwa = 51;

  auto r = lmdif1(gaussian_cminpack, &data, 9, 3, iGuess, fvec, 0.01, iwa, wa, lwa);

I get r = 1 and iGuess is not changed. Why?

Also, I've tried to use lmder1, when I should give the jacobian 9x3 matrix and also I don't get some adequate result. Can you help me please?

devernay commented 3 years ago

remove abs from fvec: it has to be differentiable everywhere. also, ftol is way too high: use sqrt(dpmpar(1)) as in the example: https://github.com/devernay/cminpack/blob/master/examples/tlmdif1c.c