Conflicting covar results using `leastsq` vs `least_squares` when chi-square is close to zero

[mdaeron]

First Time Issue Code

Yes, I read the instructions and I am sure this is a GitHub Issue.

Description

Standard errors for the best-fit parameters reported by lmfit.minimize() using method='leastsq' and scale_covar=False are inaccurate when chi-square is very close to zero (i.e. when the model fits the data almost perfectly). Below I perform a straight-line regression of Y = [0,1,2] against X = [0,1,2], and add small, increasing amounts of (pseudo-)noise f to the data. Only when the noise is large enough are the reported uncertainties equal to those expected from ordinary least squares.

By contrast, using method='least_squares' yields the correct standard errors. As far as I know, the Levenberg-Marquardt ('leastsq') and Trust Region Reflective ('least_squares') methods differ in how they approach chi-square minimum but should use the same underlying variance-covariance calculations after that. I thus suspect there is a bug in the leastsq method's code.

A Minimal, Complete, and Verifiable example

import sys
import numpy as np
from lmfit import minimize, Parameters, __version__

print(f'Using Python {sys.version}, lmfit {__version__}')

def residual(p, x, y):
    return y - (p['a'] + p['b'] * x)

def report(out):
    print(f'{f:8.0e} {out.chisqr:8.0e} {out.params["a"].value:8.2f} {out.covar[0,0]**.5:8.2f} {out.params["b"].value:8.2f} {out.covar[1,1]**.5:8.2f} {out.covar[0,1]/(out.covar[0,0]*out.covar[1,1])**.5:8.2f}')

params = Parameters()
params.add('a', value=0.5)
params.add('b', value=0.5)

for method in ['least_squares', 'leastsq']:

    print(f'\nMETHOD = {method}\n')
    print(f'{"noise":>8} {"chisq":>8} {"a":>8} {"SE(a)":>8} {"b":>8} {"SE(b)":>8} {"cor(a,b)":>8}')
    print(f'{"-"*8} {"-"*8} {"-"*8} {"-"*8} {"-"*8} {"-"*8} {"-"*8}')

    f = 0
    X = np.array([0.,1.,2.])
    Y = np.array([0.,1.,2.])
    out = minimize(residual, params, args=(X, Y), scale_covar = False, method = method)
    report(out)

    for e in np.linspace(16,2,15):
        f = 10**-e
        Y[1] = 1 + f
        out = minimize(residual, params, args=(X, Y), scale_covar = False, method = method)
        report(out)

Output:

Using Python 3.7.6 (default, Jan  8 2020, 13:42:34) 
[Clang 4.0.1 (tags/RELEASE_401/final)], lmfit 1.0.2.

METHOD = least_squares

   noise    chisq        a    SE(a)        b    SE(b) cor(a,b)
-------- -------- -------- -------- -------- -------- --------
   0e+00    2e-32    -0.00     0.91     1.00     0.71    -0.77
   1e-16    2e-32    -0.00     0.91     1.00     0.71    -0.77
   1e-15    8e-31     0.00     0.91     1.00     0.71    -0.77
   1e-14    7e-29     0.00     0.91     1.00     0.71    -0.77
   1e-13    7e-27     0.00     0.91     1.00     0.71    -0.77
   1e-12    7e-25     0.00     0.91     1.00     0.71    -0.77
   1e-11    7e-23     0.00     0.91     1.00     0.71    -0.77
   1e-10    7e-21     0.00     0.91     1.00     0.71    -0.77
   1e-09    7e-19     0.00     0.91     1.00     0.71    -0.77
   1e-08    7e-17     0.00     0.91     1.00     0.71    -0.77
   1e-07    7e-15     0.00     0.91     1.00     0.71    -0.77
   1e-06    7e-13     0.00     0.91     1.00     0.71    -0.77
   1e-05    7e-11     0.00     0.91     1.00     0.71    -0.77
   1e-04    7e-09     0.00     0.91     1.00     0.71    -0.77
   1e-03    7e-07     0.00     0.91     1.00     0.71    -0.77
   1e-02    7e-05     0.00     0.91     1.00     0.71    -0.77

METHOD = leastsq

   noise    chisq        a    SE(a)        b    SE(b) cor(a,b)
-------- -------- -------- -------- -------- -------- --------
   0e+00    2e-24     0.00     1.00     1.00     0.45    -0.00
   1e-16    2e-24     0.00     1.00     1.00     0.45    -0.00
   1e-15    2e-24     0.00     1.00     1.00     0.45    -0.00
   1e-14    1e-24    -0.00     0.00     1.00     0.50    -0.45
   1e-13    2e-24    -0.00     0.00     1.00     1.00    -0.89
   1e-12    2e-24    -0.00     0.00     1.00     0.50    -0.45
   1e-11    7e-23     0.00     1.00     1.00     0.45    -0.00
   1e-10    7e-21     0.00     1.00     1.00     0.45    -0.00
   1e-09    7e-19     0.00     1.00     1.00     0.45    -0.00
   1e-08    7e-17     0.00     0.24     1.00     0.50    -0.44
   1e-07    7e-15     0.00     0.93     1.00     0.67    -0.74
   1e-06    7e-13     0.00     0.92     1.00     0.70    -0.77
   1e-05    7e-11     0.00     0.91     1.00     0.71    -0.77
   1e-04    7e-09     0.00     0.91     1.00     0.71    -0.77
   1e-03    7e-07     0.00     0.91     1.00     0.71    -0.77
   1e-02    7e-05     0.00     0.91     1.00     0.71    -0.77

Version information

Python: 3.7.6 (default, Jan  8 2020, 13:42:34) 
[Clang 4.0.1 (tags/RELEASE_401/final)]

lmfit: 1.0.2, scipy: 1.4.1, numpy: 1.18.1, asteval: 0.9.18, uncertainties: 3.1.4

Link(s)

https://stackoverflow.com/questions/66883326/inaccurate-parameter-uncertainties-from-lmfit-when-chi-square-is-close-to-zero

[newville]

@mdaeron Please explain why is this a bug in the code or close this issue and start a conversation on the mailing list, as (really, no kidding) clearly instructed.

It seems to me that since you linked to a question on the same topic on Stackoverflow, that you know that it is a question for discussion and clarification, which is exactly the situation for which we say "start a conversation on the mailing list, not here as an Issue in the libraries bug tracker". It is hard to believe that you did not understand this.

We've said this really so many times and in so many different ways that it is HUGE demotivator when people refuse (yes, I do mean that word to convey a deliberate decision) to understand this.

Sure, I suppose we could just ignore all of our recommendations for how to ask questions. This is free software written and maintained entirely by volunteers who have other jobs and interests. We have maintained and documented the library and gone out of our way to set up forums to try to help people to use these tools You sort of ignored all of that.

[mdaeron]

@newville when I asked the question on SE I was unsure if this was a bug. Since then, I discovered that least_squares provides the correct values. Because least_squares and leastsq appear (going by the lmfit documentation) to use the same approach to compute variance-covariance, the fact that the latter yields inaccurate results very reasonably suggests that this is indeed a bug, so I opened the present issue (adding the comparison between methods). If, on the other hand, this is an (undocumented) feature or if you're just not interested, please just close the issue.

At any rate, I stand by my initial bug report. YMMV, of course. Volunteers (who have other jobs and interests) reporting bugs is a feature of FOSS in my opinion.

Thanks again for lmfit, and have a nice day.