FuncFit and FindRoots

Yes, you can use FindRoots in a fitting function, but you have to be careful. The solution from FindRoots has limited precision, and may not give absolutely smooth output. Consequently, you need to use an epsilon wave. First I would run your fit function outside of a curve fit to verify that it actually computes reasonable values.

Here is a description of the epsilon wave that I wrote for another customer:

Curve fitting uses partial derivatives of your fit function with respect to the fit coefficients in order to find the gradient of the chi-square surface. It then solves a linearized estimate of the chi-square surface to find the next estimate of the solution (the minimum in the chi-square surface).

Since Igor doesn't know anything about your fit function, it (he?) must use a numerical approximation of the derivatives, which are calculated by finite differences. That is, a model value is calculated at the present estimate of the fit coefficients, then each coefficient is perturbed by a small amount and the derivative is calculated from the difference.

The epsilon wave sets the size of the perturbation used to calculate the derivative. If you don't use an epsilon wave, Igor sets epsilon to:

if your coefficient wave is single precision,
eps[i] = 1e-4*coef[i], unless coef[i] is zero, in which case eps[i] = 1e-4
if your coefficient wave is double precision,
eps[i] = 1e-10*coef[i], unless coef[i] is zero, in which case eps[i] = 1e-10

There are a couple of reasons you might need to explicitly set the epsilon. One is if your fit function is insensitive to a coefficient. That is, perturbing the coefficient makes a very small change in the model value. You can get a situation where the dependence of the model is so small that floating-point truncation results in no change in the model. You have to set epsilon to a sufficiently large value that the model actually changes when the perturbation is applied. But this most likely not your problem.

Another case where you need to set epsilon is where the model is discrete or noisy. This can happen if the model involves a table look-up or a series solution of some sort. In the case of table look-up, epsilon needs to be large enough to make sure you get two different values out of the table.

In the case of a series solution, you have to stop summing terms in the series at some point. If the truncation of the series results in less than full floating-point resolution of the series, you need to make sure epsilon is large enough that the change in the model is larger than the resolution of the series. A series might include something like a numerical solution of an ODE (using IntegrateODE). It could also involve FindRoots or Optimize, each of which gives you an approximate result, and run faster if you don't demand high precision. Now the epsilon value must be set large enough that you get model values sufficient different that they are larger than the noise introduced by FindRoots.

In the case of FindRoots, you may also improve (but maybe not entirely solve) the problem by using /T to increase the accuracy of the result. That would be at the expense of longer computation time.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

October 18, 2011 at 11:38 am - Permalink

ChrLie

Thanks for the advice! Setting epsilon values did the job!

Cheers
C

October 19, 2011 at 02:59 am - Permalink