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

问题描述

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

这是 lmfit 如何估计参数不确定性的错误/限制，还是我对 lmfit 期望值的理解的错误/限制？

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

解决方法

作为记录，当使用 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.,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,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 python

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

问题描述

解决方法

相关问答