当卡方接近于零时，lmfit 的参数不确定性不准确

Question

当卡方非常接近于零时（即当模型几乎完美地拟合数据时），我的 lmfit（当 scale_covar = False 时）报告的最佳拟合参数似乎不准确.下面我执行了 Y = [0,1,2] 对 X = [0,1,2] 的直线回归，并向数据中添加了越来越多的噪声 (f)。只有当噪声足够大时，报告的不确定性才等于普通最小二乘法预期的不确定性。

这是 lmfit 如何估计参数不确定性的 bug/limitation，还是我对 lmfit 的预期的 bug/limitation？

# python 3.7.6
# lmfit 1.0.2

import numpy as np
from lmfit import minimize, Parameters

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}')

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}')

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

f = 0
X = np.array([0.,1.,2.])
Y = np.array([0.,1.,2.])
out = minimize(residual, params, args=(X, Y), scale_covar = False)
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)
    report(out)

# OUTPUT:
# 
# 
#    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

Answer 1

郑重声明，当使用 least_squares（信任区域反射）方法而不是默认的 leastsq (Levenberg-Marquardt) 方法时，上述行为会发生变化，而是返回准确的估计值（请参阅下面的代码）。

在许多（大多数？）情况下，前一种选择在收敛方面等同于或优于后者，并且在（边缘）情况（例如所描述的情况）中更稳健地估计方差-协方差的额外好处这里。我不确定使用 leastsq 比 least_squares 有什么明显的好处。

import numpy as np
from lmfit import minimize, Parameters, __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:
# 
# 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

当卡方接近于零时，lmfit 的参数不确定性不准确

Inaccurate parameter uncertainties from lmfit when chi-square is close to zero

python

lmfit