I have made a bunch of changes to gnuplot fitting in fit.c.
Generally, they are flagged by comments containing "tsm feb06".
In many cases, old code is still there for reference but commented out.
This should probably be cleaned up before a release.
Some of the changes have user-variables to control them.
I have not yet figured out how to commit the changes,
so they aren't out there yet.
1. Error message for undefined function evaluation
The old message in call_gnuplot() just said "Undefined value during
evaluation." I added the data point number, x, y, z, and parameter
2. Check for error value equal to zero
If the user supplies errors but any of them are zero, there will be a
zero, trashing the fit. I added a check in fit_command(), which prints
point number, x, y, z, and error, then aborts back to command line.
3. Made new message when fit stops because lambda is too large
Previously, the message in regress() said that the fit had converged in
which wasn't always true. Other output unchanged.
4. Zero-change in chisquare is now "BETTER" rather than "WORSE"
Previously, if the exact minimum of chisquare was found, so the
was zero, this was called "WORSE" in marquardt(). The loop in
regress() didn't exit,
and further iterations were done, which were also "WORSE" and increased
until the upper limit on lambda or iterations was reached. This was a
time, and potentially confusing, for no benefit.
5. Final parameter cleanup fixes
The gnuplot internal variables for the fit parameters may contain values
from an iteration that made the chisquare worse rather than better, if
the loop in regress() terminates due to maximum iterations or maximum
rather than convergence. There was code at the end of regress() that
the internal variable for only the last parameter, for a different
(the last variable is altered in call_gnuplot() to calculate
but not restored there). I changed the code at the end of regress() to
_all_ internal variables, from the array that contains the parameters
the best chisquare so far. I also put a copy of the last-parameter-fix
into call_gnuplot(), where it logically belonged anyway. Without this,
user interrupted the fit and tried to plot the current function on top
of the data,
the last internal parameter was not correct.
6. Changed convergence criterion, with new user-variable FIT_LIMIT_ABS.
The new simpler criterion in regress() is absolute reduction in
for an iteration of less than epsilon*chisquare plus epsilonAbs.
for new internal variable epsilonAbs is zero, user-settable by
Default for the existing internal variable epsilon is 1.0e-4,
by FIT_LIMIT (unchanged). In most cases, the old convergence criterion
_relative_ change in chisquare of less than epsilon. But if the
was less than NEARLY_ZERO = 1.0e-30, (which could happen if fitting data
with magnitude less than 1.0e-15 with no errors), the old criterion was
_absolute_ change in chisquare less than epsilon. Any change in
would probably be less than epsilon in these cases, so the fitter would
announce convergence immediately rather than finding the minimum. Users
would probably prefer the default relative convergence criterion in
but there was no way for them to impose it. All they could do was to
FIT_LIMIT, which would only let them adjust the _absolute_ convergence
which would require them to guess what their minimum really was. With
new criterion, the default convergence criterion is always relative no
what the chisquare is, but users now have the flexibility of adding an
convergence criterion through FIT_LIMIT_ABS. The new scheme defaults to
doing the same thing as the old scheme in almost all cases, is more
in some (rare) cases, is easier to understand, and is more flexible.
The load-file newfit0.dem tries to fit a line to the file newfit0.dat.
With old gnuplot, the fit gives up when the parameters are still pretty
with the new convergence criteria the fit continues until it is right.
7. New one-line progress-report, revert by FIT_CLASSIC_PROGRESS = 1
The old report was many lines per iteration, so you could not hold many
iterations in the scroll buffer, and it was hard to see how chisquare,
parameters, and lambda were changing. It also said nothing about
where the chisquare increased (except to print an asterisk). It printed
"WSSR" with no explanation that this means weighted sum of squared
rather than the more common term "chisquare". The new report is one
per iteration, with everything in neat columns so it's easy to track
The first column is the chisquare, the second is the change in
by the convergence limit, so convergence means a value between -1 and
lambda value, and one column per parameter. It also shows iterations
the chisquare increases. Implementation is through new function
called by regress() [and now also by marquardt()]. Both console output
fit.log file are affected by the change. The user can revert to the old
progress report show_fit() by setting FIT_CLASSIC_PROGRESS = 1.
8. New final fit parameter report format, revert by FIT_CLASSIC_RESULT
The old default final report label "final sum of squares of residuals"
was misleading because if the user supplied errors, the number printed
was the chisquare (_weighted_ sum of squared residuals or WSSR). The new
report labels are "chi-squared, "degrees of freedom," and "chisq/ndf"
when the user supplies errors. If the user did not supply errors, it
"Sum of squared residuals" and "degrees of freedom", and not SSR/ndf
it is not particularly meaningful in that case. In both cases, it
an error-rescaling factor (square root of internal chisquare per degree
calculated with error=1 if errors are not supplied). A new internal
errColumn is set to 1 in fit_command() if the user supplied errors,
it is zero. If the user supplied errors, both the raw and rescaled
errors are printed. If the user did not supply errors, only the
parameter errors are printed. In both cases, the same thing goes to
file. The old printing code was moved into new function
the new printing code is in PrintResults2(). A new internal variable
controls which function is called by regress(). It defaults to 0 (new
but setting FIT_CLASSIC_RESULT = 1 restores the old results report.
The load-file newfit1.dem fits lines to the files newfit1.dat and
switching back and forth between the old and new progress and results
9. Error-rescaling control
The old default was to always rescale parameter errors by the square
of the "chisquare" per degree of freedom (unless there were zero degrees
of freedom, in which case errors were considered undefined). While
the only sensible thing to do if the user does not supply errors to the
it is not the only sensible thing to do if errors are supplied. The
errors are calculated in regress() as before, but now they are only
by the square root of chisquare per degree of freedom if new internal
rescaleErrs is true (non-zero). The rescaleErrs variable is set in
along with errColumn, but at the moment it is always set to 1 whether
supplied errors or not. This means that the internal variables for
(not default, but available through recompiling with a preprocessor
the (non-default) old-style result report always give only rescaled
errors, as before,
unless someone changes the source and recompiles. The new default
from FitResults2() checks rescaleErrs and errColumn so it can print the
(both raw and rescaled if the user supplied errors, only rescaled if
the user did
not supply errors) no matter how rescaleErrs is set.
10. Derivative-step algorithm improvement, revert by
FIT_CLASSIC_DRV_STEP = 1
The fit needs derivatives with respect to parameters of the prediction
data point. These are calculated numerically by changing the parameters
by a small
amount. The old step size was DELTA=0.001 times the parameter value
if the parameter is less than 1e-30). In some cases, this step is so
non-linearities in the function make the derivative inaccurate. This
prevents convergence, and may make the location of the minimum wrong.
cases, the step is so small that the predicted function value does not
to machine precision, so the calculated derivative is zero. This
a singular matrix error. There is an optimal step size that balances
and non-linearities (see the discussion of numerical derivatives in
Recipes). It requires an estimate of the roundoff error of the
function, and an
estimate of the second derivative. The function that we ultimately
is the chisquare. We can estimate the roundoff error from the data
the matrices that the fit accumulates anyway can be used to estimate
derivative. The calculation adapts to the data, function, and
In most cases, the step size is much less than 0.001 times the
although it is larger when it needs to be to avoid roundoff errors.
algorithm is implemented by allocating (and freeing) in marquardt() a
drvStepDA, initialized by calling new InitDrvSteps(), updated during
by calling new CalcDrvSteps()and changing calculate() to use the array.
The new algorithm is the default, but the old algorithm (re-implemented
new functions) can be restored by setting FIT_CLASSIC_DRV_STEP = 1.
The load file newfit2.dem shows a case where the new derivative step
works but the old algorithm step is too small, resulting in a singular
The load file newfit3.dem shows a case where the old derivative step
is too large so the fit does not really converge, while the new
algorithm is OK.
11. Improved treatment of parameters near zero
Previously, an initial parameter value less than NEARLY_ZERO = 1.0e-30
was silently changed to 1.0E-30. But in some cases it is natural for
parameters to be this small or smaller. Rather than silently overriding
the user's input when near zero, the new version uses any initial value
with no change. The exception is that an initial value of EXACTLY zero
is disallowed. The reason is that the numerical derivative algorithm
needs at least the order of magnitude of the parameter for
There is a message telling the user to provide a non-zero
value of the right magnitude, or the magnitude of the expected error
if the expected value is zero. The old code was removed from
the new code is CheckParms(), called at the start of regress().
The old behavior of creating a new variable with initial value 1.0
when a fit parameter doesn't yet exist is retained, but a print
was added to createdvar() advising the user that it is better to provide
a non-zero explicit value.
12. Parameter step size limit, controlled by FIT_MAX_PAR_STEP
For non-linear functions, particularly when the initial parameter
values are far
from optimum, the fit may try to change the parameters to values where
is undefined, which terminates the fit with an error (at least in the
the user gets more information about the parameter values and data
caused the undefined result!). The only thing the user could do was to
change the starting point and pray. A new function LimitParSteps()
marquardt() to scale down the steps in all parameters so the ratio of
change to its value is no larger than internal variable maxParStep.
The new user
variable FIT_MAX_PAR_STEP controls maxParStep. The default value of
is 1.5, which allows a factor of 2.5 increase per iteration, and also
for the sign
to change. Note that maxParStep values of less than 1.0 would not
parameter to change sign, which could either be useful, or dangerous.
13. Scale-independence through multiplicative lambda, revert by
The lambda parameter in the Levenberg-Marquardt algorithm reduces the
if the calculated step actually increases the chisquare. The amount of
for the different parameters depends on details of the implementation.
implementation is "multiplicative", which is dimensionless and is
parameter scale differences. The implementation in gnuplot is
makes the performance sensitive to the relative scale of parameters and
Additive lambda requires a somewhat arbitrary initial value
multiplicative lambda can be simply initialized to 0.01. The
multiplicative lambda is easy (much larger than 1 means the step sizes
much less than 1 means the step sizes are "ideal"), while it is hard to
the value of additive lambda. My preference is for therefore for
lambda, although there are certainly combinations of data, function,
point where additive lambda will work and multiplicative lambda will
implementation is through changes to marquardt() in both initialization
of lambda. The default is now multiplicative lambda initialized to
additive lambda is still available by setting FIT_CLASSIC_LAMBDA = 1.
set FIT_START_LAMBDA to change the initial lambda value in either case.
additive lambda, supplying a FIT_START_LAMBDA value overrides the
calculation, and there is no way to turn it back on without restarting
(this has always been true).
The load file newfit4.dem shows a case where a simple line fit gives
answer with additive lambda (because it thinks it is converged when it
and multiplicative lambda gives the right answer.
14. "Skeptical" chisquare for nonlinear fits, controlled by
In a nonlinear fit, the influence of a given data point on the fit
(which depends on the parameter derivatives at the data point) can be
sensitive to the values of the parameters. This can result in slow
unless the initial parameters are chosen very close to the (unknown)
An example is a function like a*exp(b*x), where the derivatives at
large x can be
orders of magnitude larger than at small x. This interferes with
because Levenberg-Marquardt refuses to take a step if the chisquare
so if it gets close to the data at high x, even tiny parameter changes
chisquare worse (unless they are balanced in a nonlinear way). These
can be reduced in some cases by adjusting the weights of data points to
out the dependence of the derivatives on the parameters. I call this
"skeptical chisquare" because it doesn't take the initial parameter
as seriously. In some cases, the convergence time goes from thousands
iterations to only a few. However, the location of the minimum of the
chisquare will be close to the minimum of the normal chisquare if the
form can reproduce the data well, but it will not necessarily be very
the function does not reproduce the data, at least over the region
Also, the errors calculated with the skeptical chisquare are not very
For these reasons, a skeptical chisquare fit should only be used to
initial values for a normal chisquare fit. There are if-statements in
that test new variable skepticalChisq, and re-scale the internal errors.
The default is 0 for normal chisquare. The user can turn it on by
FIT_SKEPTICAL = 1,
and off by FIT_SKEPTICAL = 0.
The load file newfit5.dem gives an example with a*exp(b*x) where
takes thousands of iterations with normal chisquare, even starting with
parameters close to right, but only a few with skeptical chisquare, and
final refit with normal chisquare takes only a few iterations.
15. Monte Carlo search for initial fit parameters
There are many problems where fits only converge when the initial
rather close to the solution, and it can be tedious finding such a
To make this easier, I have implemented a Monte Carlo search for the
over a user-defined range, which is used as the initial parameters for
Levenberg-Marquardt fit. The range information is input using the "via
mechanism. The "via file" parser in fit_command() was modified to
allow (but not
require) two values per line. If any lines have two values, a Monte
for starting parameter values is done. Parameters with only a single
value are held
constant for the Monte Carlo search, but still vary in the final fit.
"# FIXED" syntax for non-parameter constants was not changed. The
default Monte Carlo
search is 1000 iterations. The user can change this by setting
The default behavior is to continue with a regular fit from the best
But if FIT_MONTE_CARLO_TRIALS is negative, after (minus) that many
trials, control goes
back to the user without a fit, but with the internal parameter values
set to the
best trial. This lets the user plot the function with these values as
The random number generator is not re-seeded, so repeating the search
with the same number
of iterations on the same function and data always gives the same
The load file newfit6.dem shows sine-wave data where fits fail unless
parameter is quite close to the right value. The Monte Carlo search
does quite a good
job in 1000 iterations, which takes only seconds, and the fit converges
well from this
I compared the new fitting code results with the old version on fit.dem.
The first fit starts with both parameters initialized to zero, which is
illegal, so I changed fit.dem to initialize them to 0.001. When the
are set to revert to the old algorithm, the results are nearly
The first fit takes a few extra iterations because of maxParStep.
With the default changes to the derivative and lambda algorithms, the
results are noticably different on the hemisphere fit. I traced this to
the fact that the function has a square root, and for some of the data
points, the argument is negative. The fit only uses the real part of
the function, so it is zero for those data points. So the chisquare
actually has multiple minima, depending on how many data points give
negative square roots. The modified fit.dem now makes a plot of the
argument of the square root, and also defines a new fit function with
a Taylor expansion of the square root that is defined for all arguments.
With the modified fit function, the new and old versions give the same
result (and both converge in far fewer iterations, implying that most
of the iterations were being confused by the multiple minima).
The new version generally tends to converge in somewhat fewer
although for some of the cases in fit.dem, the new version takes more
Prof. Thomas Mattison, Dept. of Physics & Astronomy, Univ. of British
Present Address: Stanford Linear Accelerator Center
2575 Sand Hill Road, Menlo Park, CA, 94025
Building 48 (Research Office Building), Mail Station MS35
Office: ROB-231 Phone: 650-926-5342 Fax: 650-926-8522