Re: Changes to gnuplot fitting submitted to sourceforge as patch

[<br...@ph...>]

Thomas Mattison wrote:
> On 8-Mar-06, at 10:34 AM, Hans-Bernhard Bröker wrote:

[A side note: you should subscribe to gnuplot-beta if you're going to
send mail there --- if you don't, each of your submissions will sit in
limbo until I get round to approving it...]

>>> Yeah, I wondered whether that would be sufficiently portable.
>>> Is there a random generator internal to gnuplot that I should use?
>>
>> Hmm, we do have specfun.c:ranf().  That should fit the bill.
> 
> I looked at specfun.c.  There's something declared
> static double ranf(struct value *init)
> but it's not visible outside specfun.c, and I didn't know what to do
> with the struct value* .

ranf() is the real function you need.  struct value *init() is
equivalent to state argument of the thread-safe variant of drand48.
f_rand is an interface to ranf() to be called by users, i.e. the thing
they get if they "print rand()" at the command prompt.

Well.  For the moment, let's just stick to rand()/RAND_MAX for the
moment, and fix this later.  I can export a drand48() like interface
from specfun.c any time.

>> ["TM" is Thomas Mattison, HBB are my, Hans-Bernhard Broeker's, comments]

[Points not replied to should be considered as agreed upon.]

>> 4.  Zero-change in chisquare is now "BETTER" rather than "WORSE"

[...]

> I think it would be hard to construct a test case where an iteration
> produces EXACTLY zero change in chisquare, where the following iterations
> succeeded in improving the chisquare.  The only cases I can think of
> where the chisquare change would be EXACTLY zero are that the parameters
> are already at the minimum to machine precision, or the value of lambda
> is so large that the parameter steps are so small that the chisquare
> doesn't change.  In the first case, it's OK if we stop iterating.  In
> the second case, a result of WORSE will make lambda even bigger, so
> the next step will be even smaller, so it will also produce zero change
> in chisquare.

One possibility I can imagine is that a local extremum isn't a point,
but a large exactly horizontal plateau.  This would be an ill-posed
problem, of course, but since AFAIK nobody has ever found a way of
teaching users not to pose any of those, there we go.  The correct
reaction to that kind of situation, should be to find the boundaries of
that plateau and see if the slope is up or down, out there.  If the
change does that, I'm all for it.

> I have added a print statement that tells the user the
> truth when the fit has terminated because lambda hit the maximum. But
> in order for this to be useful, we need to make sure that this
> _isn't_ printed for a fit that has actually found the minimum. The
> way to do that is to make EXACTLY zero change in chisquare be BETTER
> instead of WORSE, so the loop exits then, rather than waiting for
> lambda to hit the limit.

Ah, now there's an argument even I can understand ;-).  OK, then.

>> TM: Without this, if the user interrupted the fit and tried to
>> plot the current function on top of the data, the last
>> internal parameter was not correct.
>>
>> HBB: in that case, the fit interruption handler is what needs fixing,
>> not call_gnuplot
> 
> That would be another way of fixing things, but is there any significant
> reason why calculate() should not undo the perturbation of the last
> parameter that it had just done?

The problem is that a hard interruption could, in principle, happen
anywhere inside calculate().  At least that's how I remember this being
handled in the Linux versions, where this is done directly by SIGINT.

>> 6.  Changed convergence criterion, with new user-variable FIT_LIMIT_ABS.

>> TM:
[...]

>> HBB: I have to disagree here.  Of course users could impose a
>> change: they could supply the missing error values.  Fake them by a

> But this looks like a case where it's easier to fix the code so it works
> even when somewhat abused, rather than explaining what constitutes abuse,
> why it causes the problem, and how to get the results they want.

Let's just say I'm not fully convinced that the exact usage rules of
FIT_LIMIT_ABS will be less prone to incorrect usage than error-less fits
already are in general.  Users tend to just see a knob they can turn
which will let their fits print "converged" at the end, and never look
back to find out how it works, or whether it should be used in a
particular case.  We may be giving the poor guys a gun to shoot
themselves with instead of teaching them how to fish.

>> 7.  New one-line progress-report, revert by FIT_CLASSIC_PROGRESS = 1

>> HBB: the meaning of all abbreviations is explained in the
>> documentation.  There's no excuse for user not reading the documention
>> of a tool as complex as this.  Considerable work has gone into that
>> documentation in the past.  One result of much deliberation during that
>> activity was was to call that quantity WSSR, instead of chisquare.
>> I'd rather not go through all that again.

> Ah, if they would only read documentation....

Well, making it easier to avoid reading it isn't going to help with that
problem.  So let's not make it any easier than we have to.

> Is the discussion you mention recorded anyplace?  

I don't think so.  That was a rather lengthy email exchange back in
1997, between me and Lucas Hart of "orst.edu" (no where that is).
 From a quick peek at it, Lucas' point against using "chisquare" in the
printouts was that it must be the same as the one in "chisquare
distribution" and "chisquare test".  But that's true only in some cases.

> But it's easy to change my column label to WSSR if that's really
> more universally understood, or even to SSR for no-errors fits
> and WSSR for fits with errors.

The truth is gnuplot always does weighted SSR (sometimes the weights
are just 1.0), i.e. printing WSSR is never really wrong.  It's just
confusing until people read the docs.  Which, IMHO, is actually a Good
Thing(TM).

[...]

> I considered adding line breaks, but decided against it.  They don't
> solve the readability problem for narrow consoles, compared to letting
> the text wrap.  

Line breaks with some indentation might work, though.  Maybe yet another
new parameter: FIT_LOG_WRAP_COLUMN (zero means don't wrap)?

>> 9.  Error-rescaling control

> I do advocate changing the default behavior.  It's statistically wrong to
> rescale the errors according to the chisquare, if the user provided valid
> errors.  

Well, it's statistically wrong to take seriously _anything_ the fit
prints when the chisquare/ndf is far enough away from 1 to make a
difference.  Such fits are plain any simply inacceptable.  So to some
extent, it doesn't matter at all what we do with them: any result will
be just as wrong as any other.

 From a different point-of-view, a chisq/ndf far from 1 means the data
errors don't explain the differences between data and model.  Either the
data errors are correct --- then the model is wrong.  Or the data errors
are (as they so often are) bollocks.  'fit' has no way of knowing which
is the case.  It has to favour one of them blindly, or it has to give up
right away, and just refuse to print errors at all.

> If we repeat the same experiment and fit many times, the data will
> have statistical fluctuations, the fit parameters will have statistical
> fluctuations, and the chisquare will have statistical fluctuations.  

... and the parameter errors will also have statistical fluctuations.
Which will generally be no smaller than those of chisquare itself.
So the dividing them doesn't actually increase the variation of the
reported parameter errors considerably.

>> 10.  Gnuplot-readable parameters and errors in one line in fit.log file

> The point is that I want to _append_ many fit results to a
> gnuplot-readable file, so the results from many similar fits
> could be conveniently plotted along with their errors.  The
> file from the update command doesn't seem to be appropriate for
> appending results from many fits.

Not yet.  But 'update's job is more similar to what you're doing than
that of the fit.log file.

Sticking that machine-readable data into the middle of a human-readable
data stream, from which it'll have to be extracted before it can be
used, doesn't really look like a good idea.  A new command "update
append 'myfits.dat'" or whatever would make much more sense, from a user
interface point-of-view.  For one thing, it gives the user an
opportunity to choose which fits to put into the summary data file, and
which not to.

> Can you give a reason _why_ fit.log is the wrong place?

Because its primary purpose is to be human-readable, not
machine-readable.  Because for all you know, it already contains a lot
of data the moment you start gnuplot.  fit.log is, basically, an
electronic lab notebook, not a worksheet to collect data from various
steps of a single experiment in.

>> 13.  Parameter step size limit, controlled by FIT_MAX_PAR_STEP

>> what Marquardt's lambda parameter is designed to do.  Even if it's not
>> an exact duplicate, this seems to almost beg for a fight between
>> those two mechanisms over who gets to decide how big the step sizes
>> should be...
>> Couldn't the same effect be had by a simple penalty on the chisquare
>> to steer the existing algorithm away from those regions?

> It's not a duplication of what lambda does, it addresses a different
> problem.  

But it does so in a similar manner: limiting the step size.

> On the first iteration of a nonlinear fit, particularly if
> the initial parameters are poor, the first parameter step may be
> to a region where the function is not even defined, which results in
> a failure.  

This is not particular to the first iteration.  The algorithm can come
close to the boundary of the models definition space any time.  For all
we know, the minimum itself could be exactly on the boundary.

> Lambda doesn't reliably prevent this, because it can take
> several iterations for lambda to increase far enough to limit step
> sizes by itself.

So let it take several iterations.  If necessary, help it by signalling
undefined values with a penalty on chisquare.

>> 14.  Scale-independence through multiplicative lambda,
>> revert by FIT_CLASSIC_LAMBDA = 1

>> TM: implementation is "multiplicative", which is dimensionless and is
>> insensitive to parameter scale differences.  The implementation
>> in gnuplot is "additive" which makes the performance sensitive to
>> the relative scale of parameters and errors.
>>
>> HBB: I must admit you've lost me there.  Could you pass me some
>> references about these two different lambda's?

[...]

> It's true that a sophisticated user can rescale the problem to circumvent
> the problems with additive lambda.  

It's also true that a sophisticated fitting program can do that for him,
automatically...

 > But that's never necessary with
> multiplicative lambda.  I think multiplicative is a better default,
> because it's scale-independent and it's easier to interpret the lambda 
> value.

I'll have to ponder this for a while longer.

>> 16. Monte Carlo search for initial fit parameters

>> HBB: Nice.  So all we now miss for a complete typical fitting
>> mess-of-tools is the Nelder-Meade Simplex search algorithm.

> I thought about adding simplex, but from what I've read and heard,
> its main strength is dealing with discontinuous derivatives, which
> chisquare fits don't normally have (though it might do better on the
> original hemisphere-fit demo problem....)  

"Normally" is not something we can easily rely on.  The Simplex method 
is, of course, necessary where derivatives can't be used at all, and 
it's better than MC at finding a starting point in completely unknown 
terrain because it doesn't restrict itself to a limited parameter range. 
  It worked well for me when I was still using CERN's MINUT a lot.