Standard Error Scaling for reporting errorbars

[Erin Adkins]

Long time user of lmfit, so first thanks for making such an amazing project!

In my specific example, I am using the Minimizer class with the minimize function (least squares) in my fitting code, which successfully returns the parameters with the reported standard error.  How does the reported standard error scale (or depend) on the number of floated parameters and/or the number of points in the residual vector being minimized?

The motivation behind the question is that I have been using the reported standard error for the parameter as the error bars on that value.  Talking with colleagues we got into a discussion about the 1/sqrt(N) difference between the standard deviation and standard error.  The question would be for lmfit implementation of standard error what would the N value be.  

I appreciate any and all help!

[Matt Newville]

HI Erin,

Background:

-----------------

In a Statistics 101 for Physical Scientists, you might get introduced to a data fitting problem with a definition of chi-square as

    chi_square = Sum_i={1, N}  { (data[i] - model[i, parameters]) / epsilon[i] }**2

that is there are N data points, M parameters (with M<N) and with N data points, and with epsilon being the estimate of the noise/variance in the data.

In the ideal case, a "Good Fit" will have chi_square  ~(N-M), the degrees of freedom in the fit.  So, "reduced chi-square" is 

   reduced_chi_square  = chi_square / (N-M)

which will mean that a good fit will have reduced_chi_square = 1, in the ideal case. And error bars on parameters are typically defined such that they increase chi_square by 1.  

Note that chi_square is used for two different purposes:

     a) assessing the quality of the fit

     b) estimating the error bars on the parameters.

In the field I work in there are decades-long debates about many of the terms in the definition of chi_square:

     1.  what is N?  No, really this is subtle and hard.  What constitutes an independent measurement? If you think every sample is obviously independent,  you are probably wrong.

     2. how do you accurately estimate epsilon?  Basically, no one analyzes data for which simple estimates of uncertainties dominate.   You can always make a better measurement until the point where you would spend more time estimating epsilon (which most sane people do not do) or the subtle assumption in your model start to fail.   So easily estimated uncertainties for the data could be off by a factor 2 or 10 or 50.

As far as I know,  in the physical sciences outside of "rare event particle physics", it is really rare to have data for which chi_square is easily determined and anywhere close to N-M.  Being off by factors of 100 is not unusual.  So, many people use chi_square (and related stats) to assess the quality of the fit.  For error bars, they basically say:

    "OK, let's assume this fit is 'Good'.  If I scale the estimate of epsilon should be scaled by sqrt(reduced_chi_square).   

     That would make reduced_chi_square=1, and then the error bars will be those that increase chi_square by 1"

Another way to think about that is to say the error bars are those that increase chi_square by reduced_chi_square.

What happens in lmfit

------------------------------

By default (when `scale_covar=True`) lmfit reports "chi_square" and "reduced_chi_square" according to the epsilon used (and many, many people just use epsilon=1).  The parameter uncertainties use the rescaled covariance matrix as suggested above.   This effectively asserts that epsilon should have been set such that so that reduced_chi_square = 1, and then the reported uncertainties on the parameters are the "1-sigma" uncertainties or standard error.  

For test cases where the uncertainties are actually known, this has been shown to be reliable.  Perhaps surprisingly, the simple and fast estimates of parameter uncertainties from `leastsq` are usually quite close to the values determined by more-or-less brute force with conf_interval().   You have to go out of your way or have very asymmetrically and highly correlated parameters for the uncertainties from `leastsq` to be wrong enough to be misleading. 

 

Hope that helps,

--Matt 

[Erin Adkins]

Hi Matt, 

Thanks so much for the detailed response!  The Statistics 101 primer definitely provided me with the language to discuss this nuanced issue more intelligently going forward.  

Your point about what is N and what constitutes an independent measurement was exactly the point that we are attempting to explore.  In our data the number of points in the fit may drastically increase, but for the most part that is not acting as a new independent measure that should lower parameter uncertainty.  This  gets to the point that I want to clarify from your response.  

In the what happens in lmfit portion you mention that the when the scale_covar = True (it is in our code and currently epsilon is = 1), the parameter uncertainties use the rescaled covariance matrix as suggested above.  Can you provide some more description on what that means? And specifically what the difference is between the scale_covar = True and scale_covar=False cases in terms of the reported parameter uncertainty?

Thanks again for your help!

Erin

[Matt Newville]

Hi Erin, 

On Tue, Sep 17, 2019 at 7:16 AM Erin Adkins <eadki...@gmail.com> wrote:

Hi Matt, 

Thanks so much for the detailed response!  The Statistics 101 primer definitely provided me with the language to discuss this nuanced issue more intelligently going forward.  

Your point about what is N and what constitutes an independent measurement was exactly the point that we are attempting to explore.  In our data the number of points in the fit may drastically increase, but for the most part that is not acting as a new independent measure that should lower parameter uncertainty.  This  gets to the point that I want to clarify from your response.  

In the what happens in lmfit portion you mention that the when the scale_covar = True (it is in our code and currently epsilon is = 1), the parameter uncertainties use the rescaled covariance matrix as suggested above.  Can you provide some more description on what that means? And specifically what the difference is between the scale_covar = True and scale_covar=False cases in terms of the reported parameter uncertainty?

With scale_covar = True, the parameter uncertainties are those that would increase chi_square by reduced chi_square, as if you were not highly confident in your data uncertainties (epsilon), and wanted to assert that chi_square really should equal N-M.

With scale_covar = False, the parameter uncertainties increase chi_square by 1, as if you were confident that the data uncertainties were correct.

--Matt