Re: Changes to gnuplot fitting submitted to sourceforge as patch

[Thomas M. <mat...@ph...>]

On 10-Mar-06, at 8:52 AM, Hans-Bernhard Br=F6ker wrote:

>> On 8-Mar-06, at 10:34 AM, Hans-Bernhard Br=F6ker 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...]

I did subscribe, so I don't understand.  But I have several email=20
accounts,
and some of them have aliases for the mail servers, so gnuplot-beta
may not have figured out who I am.

>
>>> 4.  Zero-change in chisquare is now "BETTER" rather than "WORSE"
>
> [...]
>> 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.

For an exactly horizontal region, the gradients of chisquare are all=20
zero,
so the calculated step will be zero.  Continuing the iteration in this
case will just result in multiple zero-length steps, until lambda maxes=20=

out.



<change: restore internal parameters from best iteration at end of=20
regress(),
  and also always restore last internal parameter inside calculate()>

>
>>> 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=20
>> 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.

It looks like the sigint handler just sets a software flag and returns,
and regress() checks the flag to know when to run fit_interrupt().

I looked, and fit_interrupt() already does restore the internal=20
variables
from what I think is the best iteration before running any script.

So I guess there was already a workaround for the fact that calculate()
left the last parameter changed, which also restored the internal=20
variables
set to the best iteration, not necessarily the one that was interrupted.

But it was true that if regress() stopped without really converging,
that the internal parameters were not left set to the best iteration
(except for the last one!)

I still think it is better if calculate() doesn't have
mysterious side-effects, even if they are fixed elsewhere.




>>> 6.  Changed convergence criterion, with new user-variable=20
>>> FIT_LIMIT_ABS.
>
> 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=20
> 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.

My perspective is that the default behavior will now be to use a fully
scale-independent relative convergence criterion, which is what I think
naive users expect, and not what they got before.  If anyone has a real
need for an absolute convergence criterion, and managed to exploit the
odd behavior of the old convergence criterion to get one, now they can
get the same thing in a straightforward way.


>>> 7.  New one-line progress-report, revert by FIT_CLASSIC_PROGRESS =3D =
1
>
>> 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=20
> another
> new parameter: FIT_LOG_WRAP_COLUMN (zero means don't wrap)?

When it wraps, the columns typically don't line up, so usually the only
thing that a line-wrap parameter would improve is that a single column
would not get broken across lines.  At that level, I'm not sure it's
worth doing.  Even with a column split between lines it's easier to use
the new  progress report to check progress than it was to use the old
format, and the old format is still available.  I think most people=20
don't
even look at the progress reports, so it's not that big an issue.



>>> 9.  Error-rescaling control
>
>> I do advocate changing the default behavior.  It's statistically=20
>> wrong to
>> rescale the errors according to the chisquare, if the user provided=20=

>> 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.

For cases where the chisquare is close to 1/DOF, it's wrong to rescale
the errors.  I do agree that when the chisquare is grossly out of whack,
it doesn't matter much what we do.

> =46rom 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=20=

> the
> data errors are correct --- then the model is wrong.  Or the data=20
> 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=20=

> up
> right away, and just refuse to print errors at all.

My solution, which you haven't complained about, is to print out
both raw and rescaled errors when the user provides data errors
in the new result format.

The remaining disagreement is what to do in the old result format,
and what about the internal variables containing the errors.

It sounds like providing another variable to control the behavior
is the appropriate solution.  So I'll provide one in the next version.


>> If we repeat the same experiment and fit many times, the data will
>> have statistical fluctuations, the fit parameters will have=20
>> 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.

I agree that for a normal chisquare behavior we are only talking about
fluctuations of the errors by of order a factor of 2, not a factor of 10
or more.  But the fit _errors_ for repeated experiments should actually
not fluctuate at all, for fixed data errors.  Only the fit parameter
_values_ should fluctuate.


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

>>> 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=20
> 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.

The simplest solution sounds like my original proposal of a new file
with the summary lines always appended, with enough comments between
them for human editing if required.  If we don't need a way to control
whether or not the present summary goes into fit.log, we don't need a
way to control whether the one-line summary goes to the new file.


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

>> 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.

Maybe I can figure out a way to have marquardt() return WORSE
if any function evaluation is undefined.  That would let the
lambda-adjustment mechanism do part of what I'm trying to accomplish
(recover from parameter guesses that cause undefined function values).

There is still another goal, which is to avoid long jumps
that might find the wrong minimum or be slow to recover from,
even if they don't cause an undefined function evaluation.
A separate explicit parameter step limit is the only way to get this.

>
>>> 16. Monte Carlo search for initial fit parameters
>>

>> 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=20=

> is, of course, necessary where derivatives can't be used at all, and=20=

> it's better than MC at finding a starting point in completely unknown=20=

> terrain because it doesn't restrict itself to a limited parameter=20
> range.  It worked well for me when I was still using CERN's MINUT a=20
> lot.

The goal of _this_ Monte Carlo is to find a starting point in a defined=20=

range,
for cases where this is hard manually, like the frequency-fit demo in=20
my patch.
Simplex is more like Levenberg-Marquardt in the sense of still needed a=20=

good
starting point.



                                         Cheers

Prof. Thomas Mattison, Dept. of Physics & Astronomy, Univ. of British=20
Columbia

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