gnuplot-beta

Re: Re: interpretation of fit errors

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

Hi

Consider the case where the "function" to be fit is just a constant: we think the the y-value of all data points should be the same, independent of x-value.  Let's call this constant C.  Surely we want gnuplot to be able to do this "fit" and give a meaningful error on C.

Then consider the special case where there is only a single data point: an x-value, a y-value, and a y-error.  The fit result should C=y, and the error should be sigma-C = sigma-y.  The fit should be "perfect" with a chisquare of zero.  The chisquare divided by the number of measurements is still zero.  

Internally, gnuplot's fit algorithm knows a number equivalent to sigma-C = sigma-y, but the reported error is that multiplied by chisquare/measurement, so the result is zero.

It is totally unreasonable to claim that "fitting" the single data point gives C with infinite precision, when the data point has a finite error.

One can make a similar argument about the case of fitting a line (two parameters, slope and intercept) to two x-y-yerror data points.  Gnuplot's fit will go exactly through both data points, with zero chisquare per measurement, so the reported errors on the slope and intercept will both be zero.  But clearly there is a range of slopes and range of intercepts that are "within one sigma" of both data points.  And internally gnuplot's fit algorithm knows those errors on slope and intercept, before it multiplies them by zero.

Every other fitting program that I have ever used, which takes into account y-errors, just reports the "internal" error and doesn't divide it by chisquare/measurement the way gnuplot does.

If the y-errors are known and are "correct", the points will scatter around the fit curve such that the residual (measurement minus curve) divided by the error has a gaussian distribution with mean of zero and sigma of 1.0.  If the "real" y-errors are twice as large as the "assumed" y-errors, the gaussian will have sigma of about 2.0, and the chisquare/measurement will be about 4.0.


However, the user is not required to provide y-errors to gnuplot.  In that case, the errors of the fit parameters are more problematic.  Gnuplot does what every other no-y-errors fitting program does: treat all the y-errors as being 1.0.  The "internal" error calculation is still done, but one does not expect it to really describe the errors of the parameters.  And one does not expect the "chisquare" per measurement to be approximately 1.0.

But we can use the "chisquare" per measurement to estimate the y-error, assuming that it's the same for all measurements.  The result is that y-error is just sqrt(chisquare/measurement).  

If we then repeated the fit using those y-errors, we do expect that the "internal" parameter errors will describe our knowledge of the parameters.

In fact, it's not necessary to actually repeat the fit.  The parameter values would come out exactly the same, and the parameter errors would just be multiplied by sqrt(chisquare/measurement).

This is in fact what gnuplot's fit algorithm does.  And for the case where no user errors are provided, I think it's (almost) the most sensible thing to do.  



So what should be done?  If it were me, I'd report both the "raw" and the "rescaled" parameter errors on the console.  The documentation should state that if the user has provided errors, the "raw" errors are the most meaningful, and if the user has not provided errors, the "rescaled" errors are the most meaningful.  This requires essentially no algorithm changes, just some print statement changes.

Perhaps at the end of the fit, one could print a recommendation to the user on which error to use, based on whether or not he provided errors, and on the chisquare per measurement.  If no errors were provided, only the rescaled errors are meaningful.  If errors were provided and the chisquare per measurement is reasonable (less than 2 perhaps), the raw errors are meaningful.  If errors were provided and the chisquare per measurement is not reasonable, then the raw errors are questionable and the rescaled errors might be more meaningful.


I would make one small additional change: chisquare should be divided by "degrees of freedom" == (measurements minus parameters) not just measurements.  Do this both for reporting STDFIT and rescaling of errors.  Dividing by degrees of freedom instead of measurements is totally standard statistical practice.

In the one-parameter fit to one data point and line-fit to two data points examples I give above, measurements-parameters is zero.  Ideally the chisquare is also zero, so ideally the rescaling factor is 0/0=NAN, which gives the "correct" result that the rescaled errors are meaningless (in these cases).  In practice the chisquare won't be exactly zero, so the rescaling factor will be INF, still telling the user that the rescaled errors are meaningless.  This is much more sensible than what gnuplot does now, which is report the errors are zero.

If we were fitting a line to 3 points instead of 2, with errors supplied, we expect the average raw chisquare to be 1, not 3.  If errors are not supplied, we expect the raw "chisquare" not to be 3*yerror^2, but only 1*yerror^2.  Dividing by (measurements - parameters) gives a significantly different (and more correct) estimate of the yerrors, and thus of the parameter errors.





                                      Cheers

Prof. Thomas Mattison           Hennings 276
University of British Columbia  Dept. of Physics and Astronomy
6224 Agricultural Road          Vancouver  BC        V6T 1Z1  CANADA
mat...@ph...         phone: 604-822-9690  fax:604-822-5324

Re: interpretation of fit errors

[Hans-Bernhard B. <HBB...@t-...>]

On 06.05.2011 19:48, Thomas Mattison wrote:

> Consider the case where the "function" to be fit is just a constant:
> we think the the y-value of all data points should be the same,
> independent of x-value.  Let's call this constant C.  Surely we want
> gnuplot to be able to do this "fit" and give a meaningful error on
> C.

Yes.

> Then consider the special case where there is only a single data
> point: an x-value, a y-value, and a y-error.  The fit result should
> C=y, and the error should be sigma-C = sigma-y.

It shouldn't.  The problem here is that by going down to a single data 
point, you've constructed a problem with zero degrees of
freedom.  That's not a fit --- it's just an equation to be evaluated.

Using fit for that is tool abuse.

> The fit should be "perfect" with a chisquare of zero.  The chisquare
> divided by the number of measurements is still zero.

But chisq is never divided by the number of measurements.  It's divided 
by the number of degrees of freedom, i.e. (number of measurements) - 
(number of parameters).  In the case at hand that's zero, STDFIT would 
be zero divided by zero.

> Internally, gnuplot's fit algorithm knows a number equivalent to
> sigma-C = sigma-y,

No, it doesn't, because it never even starts to work in that case:

gnuplot> fit a '-' u 1:2:3 via a
input data ('e' ends) > 0 10 2
input data ('e' ends) > e
          Read 1 points
          No data to fit

> but the reported error is that multiplied by chisquare/measurement,
> so the result is zero.

No it's not.  See above.


> It is totally unreasonable to claim that "fitting" the single data
> point gives C with infinite precision, when the data point has a
> finite error.

It's totally unreasonable to use fit on a single data point, period.

> One can make a similar argument about the case of fitting a line (two
> parameters, slope and intercept) to two x-y-yerror data points.

Same problem.  Zero degrees of freedom means 'fit' is the wrong tool for 
the job.

> Gnuplot's fit will go exactly through both data points, with zero
> chisquare per measurement, so the reported errors on the slope and
> intercept will both be zero.

No they won't.  Because gnuplot won't report _any_ errors in that case.

> But clearly there is a range of slopes
> and range of intercepts that are "within one sigma" of both data
> points.  And internally gnuplot's fit algorithm knows those errors on
> slope and intercept, before it multiplies them by zero.

No, it doesn't.

> I would make one small additional change: chisquare should be divided
> by "degrees of freedom" == (measurements minus parameters) not just
> measurements.  Do this both for reporting STDFIT and rescaling of
> errors.

What made you believe that's not already what's happening?

> user that the rescaled errors are meaningless.  This is much more
> sensible than what gnuplot does now, which is report the errors are
> zero.

Your view of the world might profit from an actual experiment.

Re: interpretation of fit errors

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

I stand corrected about two things:

1.  Gnuplot fits have always done the "right" thing and divided chisquare by degrees of freedom rather than measurements.  Sorry, I haven't looked at what gets printed out in a while; I mostly use my own fitting programs these days.

2.  Gnuplot will indeed refuse to fit a constant function to a single data point, or a line to two points.  Yes, I hadn't done the experiment.  Sorry.

However I regard item 2 as a bug not a feature.  It is perfectly well defined to ask the question, what is the chisquare for agreement between data and model, as a function of the parameter values, for these cases.  There is a perfectly well defined answer to that question.  There is a perfectly well defined point of minimum chisquare, which gives the best fit parameters.  And the curvature of the chisquare curve or surface gives the fit errors.  Nothing anomalous happens in the case where the number of data points equals the number of parameters.

I agree that in those two particular cases, it's not necessary to "do a fit" to find the parameters; it can easily be done by hand.  But mathematically the fit is still well-defined.

I hope we can all agree that when fitting a line to two data points that have errors, the slope and intercept DO have errors.  But calculating them by hand is not as trivial as calculating the slope and intercept by hand.  The easiest way to do that in practice is to run a fitting program.  The "raw" errors from the fit give the answer.

But gnuplot refuses to do that.  And the only reason I can see for the refusal is that the "rescaled" errors wouldn't make sense in that case (while the raw errors would make sense, and are useful).


Let's say I have several equal-size data sets, described by the same model, with the same amount of random noise in each.  If I fit each data set, the parameters will vary (within the errors), but they all have equal amounts of information, so they SHOULD all give the same fit errors.  The "raw" fit errors WILL be the (neglecting effects from nonlinearities which are usually small).  But the "rescaled" fit errors will vary randomly, depending on the chisquare of the fit.  For a case where someone has worked hard to make the input errors meaningful, it's wrong for the parameter error output to depend on randomness in the data.

Is it really so much to ask to have both raw and rescaled errors printed?



                                      Cheers

Prof. Thomas Mattison           Hennings 276
University of British Columbia  Dept. of Physics and Astronomy
6224 Agricultural Road          Vancouver  BC        V6T 1Z1  CANADA
mat...@ph...         phone: 604-822-9690  fax:604-822-5324

Re: interpretation of fit errors

[Daniel J S. <dan...@ie...>]

On 05/06/2011 04:10 PM, Thomas Mattison wrote:
> I stand corrected about two things:
>
> 1.  Gnuplot fits have always done the "right" thing and divided
> chisquare by degrees of freedom rather than measurements.  Sorry, I
> haven't looked at what gets printed out in a while; I mostly use my
> own fitting programs these days.
>
> 2.  Gnuplot will indeed refuse to fit a constant function to a single
> data point, or a line to two points.  Yes, I hadn't done the
> experiment.  Sorry.
>
> However I regard item 2 as a bug not a feature.  It is perfectly well
> defined to ask the question, what is the chisquare for agreement
> between data and model, as a function of the parameter values, for
> these cases.  There is a perfectly well defined answer to that
> question.  There is a perfectly well defined point of minimum
> chisquare, which gives the best fit parameters.  And the curvature of
> the chisquare curve or surface gives the fit errors.  Nothing
> anomalous happens in the case where the number of data points equals
> the number of parameters.

It seems like it should.  From H.B.B.'s last post, if there are only two 
points and the curve to fit is a line (two parameters, first order 
equation**) then the fit is exact in the sense that the residuals are 
zero...if no other information is given about the data.  But that 
doesn't say anything about the measurement "error" in the fit.  (I sort 
of prefer the word "uncertainty", but that's just me.)  We can only 
glean information about the uncertainty of the fit (i.e., the statistics 
of the original data) if there are additional data points beyond the 
number of parameters.

** Polynomial order of the fitting curve is important.  For example, 
exp() has infinite polynomial order, which has ramifications on uniqueness.


> I agree that in those two particular cases, it's not necessary to "do
> a fit" to find the parameters; it can easily be done by hand.  But
> mathematically the fit is still well-defined.

OK, so I'm beginning to see there is this vital piece of information 
about the quality of fit, reflecting the underlying statistics or the 
uncertainty in the data (i.e., how noisy it is).  That's what is meant 
by error, right?

In this discussion

http://www.erikburd.org/projects/pitfall/evaluate.html

is a column for "standard error".  Is that what is being referred to, 
the standard error typically used in statistical analysis?  (There are 
the underlying actual measurement errors, but we can't know what those 
are because there is uncertainty in the fit.)


> I hope we can all agree that when fitting a line to two data points
> that have errors, the slope and intercept DO have errors.

Yes, but I don't see how one can estimate the size of that error without 
additional data beyond two data points, unless some information is known 
apriori.


> Is it really so much to ask to have both raw and rescaled errors
> printed?

I'm fine with that idea.  But I'm surprised that no well defined, 
commonly understood terminology has developed for what seems like an 
important idea in the numerical analysis community.  Bastian points out 
that:

* Minuit does not scale
* Origin offers both options, the default being version dependent
* SAS and Mathematica default to scaling of errors

That inconsistency is the problem.

Dan

Re: interpretation of fit errors

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

On 2011-05-06, at 3:56 PM, Daniel J Sebald wrote:

> On 05/06/2011 04:10 PM, Thomas Mattison wrote:
>> 
>> 2.  Gnuplot will indeed refuse to fit a constant function to a single
>> data point, or a line to two points.  Yes, I hadn't done the
>> experiment.  Sorry.
>> 
>> However I regard item 2 as a bug not a feature.  It is perfectly well
>> defined to ask the question, what is the chisquare for agreement
>> between data and model, as a function of the parameter values, for
>> these cases.  There is a perfectly well defined answer to that
>> question.  There is a perfectly well defined point of minimum
>> chisquare, which gives the best fit parameters.  And the curvature of
>> the chisquare curve or surface gives the fit errors.  Nothing
>> anomalous happens in the case where the number of data points equals
>> the number of parameters.
> 
> It seems like it should.  

But it doesn't, from a mathematical point of view.

> From H.B.B.'s last post, if there are only two points and the curve to fit is a line (two parameters, first order equation**) then the fit is exact in the sense that the residuals are zero...if no other information is given about the data.  But that doesn't say anything about the measurement "error" in the fit.  (I sort of prefer the word "uncertainty", but that's just me.)  We can only glean information about the uncertainty of the fit (i.e., the statistics of the original data) if there are additional data points beyond the number of parameters.

You are right that there need to be more measurements than parameters in order to use the chisquare/degree-of-freedom to infer something about the errors of the fit INPUTS from the quality of the fit.

But in many cases, the user ALREADY KNOWS what the input data measurement errors are.  He knows the precision of the markings on his ruler, or how many decimal places his voltmeter has, or square-root-of-N for the bin contents of a histogram, etc [I am not discussing "systematic" errors here, like is the calibration of your voltmeter correct].

In those cases, the "raw" fit errors are the appropriate ones, and the "rescaled" errors are less appropriate.

> 
> OK, so I'm beginning to see there is this vital piece of information about the quality of fit, reflecting the underlying statistics or the uncertainty in the data (i.e., how noisy it is).  That's what is meant by error, right?

The way I would phrase it is, the noise in the input data propagates to noise in the fit parameters.  The "standard errors" of the fit parameters are the standard deviation of the input data (assumed to be known) propagated to the standard deviation of the fit parameters.

> 
>> I hope we can all agree that when fitting a line to two data points
>> that have errors, the slope and intercept DO have errors.
> 
> Yes, but I don't see how one can estimate the size of that error without additional data beyond two data points, unless some information is known apriori.
> 

Sure, but when the user has supplied meaningful y-errors on his data, we have precisely the a-priori informaton that your intuition is telling you is needed.


> But I'm surprised that no well defined, commonly understood terminology has developed for what seems like an important idea in the numerical analysis community.  Bastian points out that:
> 
> * Minuit does not scale
> * Origin offers both options, the default being version dependent
> * SAS and Mathematica default to scaling of errors
> 
> That inconsistency is the problem.

To me, that looks like most programs give you a choice.  Mathematica may default to scaling, but it looks like it's optional.  I don't have much knowledge about SAS, but I would expect a big package like that to have an option buried somewhere.

Minuit comes from the physics community, nuclear and particle physics in particular.  In physics, we tend to have pretty good a-priori estimates for the errors of the measurements going into a fit.  We use the chisquare/DOF as a measure of the (statistical) fit quality.  Physicists would seldom rescale fit parameter errors by sqrt(chisq/DOF).  And since Minuit is normally used in an expert-programming context instead of an ignorant-GUI-user context, the authors would expect the user to be able multiply by sqrt(chisq/DOF) himself.

In bioscience, and even more so in medical and social science, there's often next to nothing known about the errors of the input measurements a-priori.  In such cases, there's not much the user can do except set all the errors equal (or not tell gnuplot any errors, which will have the same result).  Then the "raw" errors from the fit aren't very meaningful, but the "rescaled" errors have at least some utility.  So defaulting to rescaled makes some sence.

My guess would be that there are more gnuplot fit users who don't know much about their input data errors than those who know a great deal about their input data errors.  So if I were forced to make a choice of only one error to present, I'd probably present the rescaled error.  It will be the "right thing" for users who don't give errors or give uniform and arbitrary errors.  And for the people who supply accurate errors, if the model really reflects the data, then chisq/DOF will be close to 1, and the "rescaled" errors will not be very different, very often, from the "raw" errors.


But why not give both?





                                      Cheers

Prof. Thomas Mattison           Hennings 276
University of British Columbia  Dept. of Physics and Astronomy
6224 Agricultural Road          Vancouver  BC        V6T 1Z1  CANADA
mat...@ph...         phone: 604-822-9690  fax:604-822-5324