gnuplot-beta

Changes to gnuplot fitting submitted to sourceforge as patch

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

Hi

My changes to gnuplot fitting are now available as a patch at 
sourceforge.

Changes are limited to src/fit.c, demo/fit.dem, and a new directory
demo/newfit, containing details in README.newfit, and many demos.

At present, the patch is constructed more for an audience of developers
to test, rather than for users to install, since the idea is for people
to evaluate it for acceptance into the distribution, rather than for
users to patch a released version.

Please run the demos in demo/newfit to see the effect of the changes.

Re: Changes to gnuplot fitting submitted to sourceforge as patch

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

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

>
> drand48() cannot be used like that.  It's quite a non-portable=20
> function.  Of the compilers I tried so far, only the most dedicated=20
> impersonators of UNIX (DJGPP and real Cygwin) have it.

On 9-Mar-06, at 5:56 AM, Hans-Bernhard Br=F6ker wrote:

> Thomas Mattison wrote:
>
>> 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* .
There's also in specfun.h something declared
void f_rand __PROTO((union argument *x))
but returning void doesn't help me much.

> ["TM" is Thomas Mattison, HBB are my, Hans-Bernhard Broeker's,=20
> comments]
>
> 4.  Zero-change in chisquare is now "BETTER" rather than "WORSE"
>
> TM: Previously, if the exact minimum of chisquare was found, so the
> chisquare change 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
> lambda until the upper limit on lambda or iterations was
> reached.  This was a waste of time, and potentially confusing,
> for no benefit.
>
> HBB: I'm not entirely sure about this one.  There could be other
> reasons for a failure to find an improvement besides having found a
> well-defined optimum.  An extended flat minimum, or a pair of local
> minima, say.  Giving up the Marquardt process immediately when the
> Newton iteration admits defeat seems to defeat the purpose of doing
> the Marquardt.

I think it would be hard to construct a test case where an iteration
produces EXACTLY zero change in chisquare, where the following=20
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.

But it is fairly easy to construct a case where the present code does
something that would likely confuse a user. As the iteration approaches
the minimum, the calculated parameter changes _should_ get smaller and
smaller, and eventually get so small that the chisquare does not change
to machine precision.  Gnuplot's default convergence limit of 1e-5
relative change in chisquare is pretty loose, so convergence is=20
typically
declared before this happens.  But if you reduce FIT_LIMIT to much
smaller values, you will find that in most cases, the fit gets to a=20
point
where further iterations no longer change the chisquare.  If this is=20
defined
as WORSE, then the iteration would continue forever.

Obviously there would be many complaints if this actually happened.  The
reason it doesn't happen is that there is _another_ way to exit that=20
loop.
When the chisquare is WORSE, then lambda is increased, and the loop=20
exits
on maximum lambda value.  The existing code does NOT print a message=20
that
the maximum lambda has been exceeded.  What it prints is that the fit=20
converged.

I know from experience teaching with gnuplot that lots of users make=20
mistakes
that cause "convergence" to nonsense.  Ideally, if gnuplot knows that=20
the fit
hasn't converged properly, it should tell the user that.  Having lambda=20=

hit
the maximum should be a sign that the fit has not converged.  I have=20
added a
print statement that tells the user the truth when the fit has=20
terminated
because lambda hit the maximum.  But in order for this to be useful, we=20=

need
to make sure that this _isn't_ printed for a fit that has actually=20
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.


> 5.  Final parameter cleanup fixes
>
> TM: 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 lambda rather than convergence.  There was
> code at the end of regress() that restores the internal variable
> for only the last parameter (for a different reason: the last
> variable is altered in call_gnuplot() to calculate derivatives,
> but not restored there).  I changed the code at the end of
> regress() to restore _all_ internal variables, from the array
> that contains the parameters that gave the best chisquare so
> far.
>
> HBB: this feels risky, at first sight.  regress() keeps more state =
than
> just the parameters --- restoring the parameters could leave them out=20=

> of
> synch with the rest, e.g. the correlation matrix and chisquare.

I looked again and I think it's OK.  The parameters array a[] and the
derivatives matrix C[] that is used for the error and correlation=20
calculation
in regress() are calculated inside marquardt().  They are always in=20
synch
with each other, and what regress() has is always from the best=20
chisquare
that marquardt() has seen, and the chisq variable is that value.


> TM: I also put a copy of the last-parameter-only fix code
> into call_gnuplot(), where it logically belonged anyway.
>
> HBB: I don't think that's correct.  call_gnuplot is not supposed to=20
> know how
> fit's numeric derivation method works.  It should just plug in the
> parameter values and evaluate.

You're right, because I made a mistake in my note, not the actual code.
I actually put the last-parameter-only fix code into calculate(), not
call_gnuplot().  calculate() is the function that perturbs the internal
variables for the derivative calculation, not call_gnuplot.  The change
is just to undo the perturbation for the last parameter, just like
is done for all the other parameters.

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


>
> 6.  Changed convergence criterion, with new user-variable=20
> FIT_LIMIT_ABS.
>
> TM:
> [...]
> less than epsilon.  But if the chisquare was less than
> NEARLY_ZERO =3D 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 chisquare 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
> these cases,  but there was no way for them to impose it.
>
> HBB: I have to disagree here.  Of course users could impose a
> change: they could supply the missing error values.  Fake them by a
> constant, if they must.  I know I tend to be a bit harsh on users
> in this regard:  I think there is a threshold of how easy it should
> be made to abuse a tool.  The above change is still OK, but IMHO
> it's getting a bit close to that threshold.


You're right that there is _something_ that they could do.  They could
either supply sensible errors so they have a real chisquare which would
not be 1e-30.  Or they could re-scale their data even without errors
so the "chisquare" came out in the range that works.

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=20
abuse,
why it causes the problem, and how to get the results they want.
We advertise gnuplot fits as working even if errors are not supplied,
and this is a case where they don't work properly.  It's easily fixed,
and the change just simplifies both the code and the explanation of
what the code does.  So why not?


>
> 7.  New one-line progress-report, revert by FIT_CLASSIC_PROGRESS =3D 1
>
> TM: (except to print an asterisk).  It printed  "WSSR" with no
> explanation that this means weighted sum of squared residuals,
> rather than the more common term "chisquare".
>
> 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....

Is the discussion you mention recorded anyplace?  I'd like to see the
arguments.  In my experience (undergrad physics as student and
professor, grad school and beyond in experimental particle physics,
sound's a lot like you actually), chisquare is by far the more
common name.

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.


> TM: The new report
> is one line per iteration, with everything in neat columns so
> it's easy to track the progress.
>
> HBB: ... but only for fits with less sufficiently few parameters to
> all fit in one line.  Please try to keep in mind that not everybody
> uses text consoles 200 or more characters wide...
>
> TM: The user can
> revert to the old progress report show_fit() by setting
> FIT_CLASSIC_PROGRESS =3D 1.
>
> HBB: Good --- because otherwise, all those users out there who have
> scripts parsing fit.log and/or console output of gnuplot would cause
> a mighty ruckus.

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.  And they make the output less useful for people who do
use wide consoles and/or capture the output for display in some other
program that can scroll horizontally.  And remember that the old multi-
line output has a different problem:  there is a limit to how many
iterations will be visible on a terminal or scroll-buffer.  To me,
the old progress report was very hard to use as a progress report,
and I think most people only looked at the final results.  The new
progress report is more likely to be useful as a progress report,
even if it's not perfect for all users, which nothing ever is.


>
> 8.  New final fit parameter report format, revert by=20
> FIT_CLASSIC_RESULT =3D 1
>
> TM: 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"
>
> HBB: see above note about nomenclature and documentation.

I'm open to debate about what to call things.  But part of my
point was that what we print means different things if the
user supplies  errors vs not supplying errors.  So we should
have different labels in the different cases.


> 9.  Error-rescaling control
>
> TM: user supplied errors or not.  This means that the internal
> variables for parameter errors (not default, but available
> through recompiling with a preprocessor option)
>
> HBB: Assuming you're referring to GP_FIT_ERRVARS there:
> that build option is supposed to be turned on by default, these days.
>
> TM: and the
> (non-default) old-style result report always give only rescaled
> errors, as before, unless someone changes the source and
> recompiles.
>
> HBB: this could be quite bad.  Changing the meaning of a=20
> machine-readable
> entity like the parameter error variables is a recipe for utter=20
> confusion.
> Please consider introducing a different set of error variables for the
> new-style errors.  There has to be a way for users'
> scripts to get at the exact same values they used to get, or at least
> to detect that the numbers now stored in the same places now mean=20
> something
> different.

I guess I wasn't clear.  My code does _not_ (presently) change the=20
values
calculated for the user-accessible variables created when GP_FIT_ERRVARS
is turned on.  The calculated errors are re-scaled using the actual
fit "chisquare" like they were before.  And if the user reverts to the
old-style final result report, they will get the same rescaled error
values as before.

I only included a comment to make it obvious how to change the source
code so the GP_FIT_ERRVARS error variables would only be rescaled by
the chisquare when the fit is done without user errors, which I think
is the right thing.  Another solution would be to create another
conditional-compile flag, but I didn't do it that way.

I don't advocate doubling the number of fit-error variables created,
I think that would be much more confusion than it's worth.

I do advocate changing the default behavior.  It's statistically wrong=20=

to
rescale the errors according to the chisquare, if the user provided=20
valid
errors.  If we repeat the same experiment and fit many times, the data=20=

will
have statistical fluctuations, the fit parameters will have statistical
fluctuations, and the chisquare will have statistical fluctuations. =20
But it's
not true that the fits that with lower chisquare are more closely=20
clustered
around the true parameter values.  They should not be assigned smaller=20=

errors.

Apparently some commercial fitting packages make this mistake also.
For instance, see http://physics.gac.edu/~huber/fitting/aapt2001.ppt
But it's still a mistake.

For the new default final result report, both raw and rescaled errors
are printed, and the code correctly prints both errors labelled properly
whether the rescaling variable is set to old-style GP_FIT_ERRVARS=20
values,
or to what I consider to be the right calculation.


> 10.  Gnuplot-readable parameters and errors in one line in fit.log =
file
>
> HBB: I'm against this.  fit.log is the wrong place for this data.
> It's much more similar to what the 'update' command already does, so =
if
> this feature is added, it should go there, not into fit.log.

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.

Can you give a reason _why_ fit.log is the wrong place?
Recall that in the exchange we had in the fall, I suggested making a
new fit-summary output file that would contain _only_ lines like this,
and your answer was

> I think it would make a lot more sense to add a couple lines to=20
> 'fit.log'
> instead of creating what would be an almost complete copy of all of=20
> its content.

So would it be best to go back to my original proposal
of a separate new file that contained only appended
one-line summaries with errors?


> 13.  Parameter step size limit, controlled by FIT_MAX_PAR_STEP
>
> TM: function LimitParSteps() called by marquardt() to scale down the
> steps in all parameters so the ratio of any parameter change to
> its value is no larger than internal variable maxParStep.
>
> HBB: It's been a while since I rewrote this Marquardt-Levenberg
> code, and I don't have my textbooks at hand,  but:
> this 'limit the step size' looks awfully like a duplicate of
> 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.  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.  Lambda doesn't reliably prevent this, because it can take
several iterations for lambda to increase far enough to limit step
sizes by itself.

For an example, try fitting my demo/newfit/newfit5.dat
with a*exp(b*x) with starting parameter values a=3Db=3D1 with old =
gnuplot.
You will get "Undefined value during function evaluation".

My modified version of gnuplot would print the extra information
that a=3D229 and b=3D7869 for the next trial, which is why the =
calculation
of the function blows up.  But with the default FIT_MAX_PAR_STEP =3D =
1.5,
it increases the value of a and b slowly enough that the function is
never undefined (although it still converges VERY slowly unless you
use FIT_SKEPTICAL =3D 1)

The two mechanisms cooperate nicely.  The new mechanism doesn't=20
interfere
with convergence once the steps are small.  It just keeps things from
running away before a reasonable value of lambda has been established
by the iteration.

The chisquare already penalizes the regions we are talking about,
the problem is that sometimes first steps are so large that the function
can't even be evaluated to calculate the chisquare!


> 14.  Scale-independence through multiplicative lambda,
> revert by FIT_CLASSIC_LAMBDA =3D 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?

I don't have a textbook or article reference, other perhaps than
Numerical Recipes (which uses multiplicative and doesn't mention
the other possibility).  Basically, you write the chisquare sum,
take the derivatives with respect to the parameters, and set them
to zero, to get N nonlinear equations in the N unknown parameters.
Then you linearize those equations, resulting in an N by N matrix
which I call the weight matrix, which is the inverse of the covariance
matrix.  The weight matrix isn't explicitly calculated in the gnuplot
implementation, but it's C times its transpose.  For a linear problem
the parameter step is the inverse of the weight matrix times the vector
of sum of residuals over errors times parameter-derivatives, and
iteration is not needed.  For a nonlinear problem, the resulting
parameter step is not optimal, and may increase the chisquare.

The Levenberg-Marquardt trick is to increase the diagonal
of the weight matrix (the diagonal elements are always positive),
so the elements of the inverse are smaller, and the parameter steps
are smaller.  If the diagonal is increased so far that the off-diagonal
elements are negligible, then the parameter steps are in a direction=20
that
is guaranteed to decrease the chisquare, but the steps will be small.

Additive lambda means add the same number lambda to each element
of the diagonal of the weight matrix.  Multiplicative lambda means
multiply each diagonal element by (1+lambda).  For multiplicative=20
lambda,
lambda=3D1000 would result in each parameter changing by 1/1000 of
the distance from its present value toward the value that would minimize
chisquare if all other parameters were held fixed for a linear problem.

If the elements on the diagonal are all about the same size, it makes
little difference whether you use additive or multiplicative lambda.
But when the elements are very different, and lambda is much bigger
than some elements and much less than some others, some parameters
are not affected by lambda and converge quickly, while others are
frozen by lambda and don't change.  So you can have false convergence.

> =46rom what I'm aware
> of the, usual approach against parameter scale differences is not
> a change to the handling of lambda, but of the parameters.  I.e.=20
> instead
> of minimizing with respect to the actual problem parameters, all=20
> parameters
> are modelled as factor*initial_value, with the factors all starting at=20=

> 1.0,
> and only the factors are minimized by the algorithm.  This makes them=20=

> all the
> same scale initially.  The only drawback is that the fit is no longer
> idempotent, i.e. re-running the fit with the converged results from an
> earlier fit can yield different results.  That's what MINUIT does, as
> far as I remember from reading its docs.

It's true that a sophisticated user can rescale the problem to=20
circumvent
the problems with additive lambda.  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=20=

value.

But it is true that some problems just seem to converge faster with
additive lambda than multiplicative lambda, so it's useful to keep both.



> 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.
> Oh, and about 100 pages of documentation explaining all this to
> gnuplot users who don't have anywhere near the experience it takes
> to avoid hurting themselves with such a mighty tool chest (whether
> by shooting themselves in the foot, or by dropping the chest on =
it...).


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....)  But the gain didn't seem
to be big enough often enough to be worth creating the new controls.

But the Monte Carlo search solves a more common problem: finding
a reasonable starting point, particularly if there can be multiple
minima.  It needs extra user input for the ranges, so I hesistated,
but realized that the parameter-file mechanism could let me add it
in a way that wouldn't burden people who didn't need/want it.

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.

Re: Changes to gnuplot fitting submitted to sourceforge as patch

[Bastian M. <bma...@we...>]

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

Personally I would prefer new options to `set fit` instead of dozens of n=
ew
FIT_xxx variables. This would be more consistent with setting other
options in gnuplot.

Attached you find an extension to gnuplot's fit which I have been using
for a while. It adds the possibility to turn of error scaling via
`set fit errorscaling` and saving of fit information (chisq, dof, WSSR)
to user variables as requested by Hans Boie (SF #1117724 [fit] access
to resulting chisquare) via `set fit fitvariables`. It also fixes
(SF #1324672 ] doc: false reference to "set fit").

These are just tiny modifications and could probably be integrated
into this large and very nice patch.

Bastian

--=20
Bastian M=E4rkisch
Physikalisches Institut, Universit=E4t Heidelberg

Re: Changes to gnuplot fitting submitted to sourceforge as patch

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

On 10-Mar-06, at 10:45 AM, Bastian Maerkisch wrote:

>
> Personally I would prefer new options to `set fit` instead of dozens 
> of new
> FIT_xxx variables. This would be more consistent with setting other
> options in gnuplot.

Controlling things through a "set fit" mechanism rather than FIT_xxx
is an appealing idea, I wish I had thought of it.  It makes controlling
fits more like controlling other things in gnuplot.  Certainly for
_new_ controls where there is no back-compatibility issue it sounds
better than more FIT_xxx clutter in the namespace.  I just copied the
code in fit.c which didn't use set, because that was easier.  I'll
look at your patch to learn how to use set and make an alternative 
patch.

Perhaps we could also have "set fit xxx" control things presently
controlled via FIT_xxx, then slowly deprecate FIT_xxx.  We'd have
to have a rule for how to resolve conflicts if both mechanisms were
trying to control things at the same time.  The simplest is just to
let the most recent change apply, and change both the set/show state
and the FIT_xxx state.  It would then be nice to issue a warning if
the value is set with one mechanism then changed using the other
mechanism, since that is more likely to be a mistake than changing
the value many times through the same mechanism.  I don't know if
dual-control would be easy or hard, I haven't looked at the 
implementations.

I'd also like to make the fit_command parser accept the same syntax
for errors that the plot command does, for further consistency 
improvement.
But that's not my highest priority.

> Attached you find an extension to gnuplot's fit which I have been using
> for a while. It adds the possibility to turn off error scaling via
> `set fit errorscaling` and saving of fit information (chisq, dof, WSSR)
> to user variables as requested by Hans Boie (SF #1117724 [fit] access
> to resulting chisquare) via `set fit fitvariables`. It also fixes
> (SF #1324672 ] doc: false reference to "set fit").
>

Having a switch like this may be the only way to get me and Hans-
Bernhard to stop going back and forth about it ;-)

And having variables for chisq, etc turned on by GP_FIT_ERRVARS
switch is probably a good idea (although then we'd have to agree on 
names.. ;-)

> These are just tiny modifications and could probably be integrated
> into this large and very nice patch.
>
> Bastian
>

It's nice to be appreciated!

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

Re: Changes to gnuplot fitting submitted to sourceforge as patch

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

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

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

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

Mailman only knows the actual mail address you subscribed as.  If you 
use something else as the FROM: field in your submissions, it'll reject 
them.  This one passed directly.

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

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

That's fixable.  E.g. if the lines are set to wrap, have them wrap 
always.  I.e. put all the non-parameter info in one line, then as many 
parameters as fit per line, with some indentation to provide some visual 
guidance.

But don't let me detain you --- if you don't implement, I might just
do it myself, once we integrated your patch.

>>>> 9.  Error-rescaling control

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

OK, then.

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

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

Now, that's a very strange statement, I think.  How could a process 
based entirely on statistically fluctuating data *avoid* fluctuation in 
some of its results?  There's a direct algebraic connection data --> 
residuals --> parameter errors.  I really don't see how the data can 
fluctuate, but the parameters not.

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

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

What I don't really like about that version is the "always".  It 
basically means the user has to re-do the entire procedure, or manually 
edit the file in the middle of a gnuplot sessions, to remove fits gone 
bad from this machine-readable log, before he forgets which fits are to 
be kept, and which not.

'update' is the existing command to, so-to-say, 'bless' a fit result as 
accepted for further usage.  That's why I think it's the right place to 
add the machine-readable session log.

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

Because the fit.log file is for humans to read, and it provides all the 
context the user can possibly need.  A full-blown copy of it in a 
machine-readable format would not serve much of a purpose that couldn't 
already be had by machine-translating the existing fit.log.   A selected 
subset of it must have the selection done inside gnuplot, not 
afterwards.  That's where update comes in handy.

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

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

But such tactics require implied knowledge about which long jumps are 
bad, and which aren't.  I don't see how it can be part of generic 
fitting program's job to second-guess the individual problem.  If the 
user already knows where the minimum is, she shouldn't be running 'fit'
to find it.

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

> The goal of _this_ Monte Carlo is to find a starting point in a
> defined range, 

I'm fully aware of that, and I'm not arguing we remove it.  FWIW, MINUIT 
also has such a method available.  But I don't think anyone I know ever 
used it on a regular basis.  A well-informed guess at startup parameters 
basically always outperformed it.

Ultimately, the statement I took from the manual for fudgit still holds: 
Non-linear least-squares fitting is an art!  It takes some learning to 
master it.